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ORBITAL PREDICTION AND DIFFERENTIAL CORRECTION 
USING VINTI’S SPHEROIDAL THEORY FOR ARTIFICIAL SATELLITES 

/[ ^3 

ABSTRACT 

The spheroidal theory developed by Vinti for determining the orbit of an 
artificial satellite of an oblate planet is presented in algorithmic form, in 
which empirically derived initial conditions are used to obtain the co-ordinate 
and velocity components of an unretarded satellite at any time. A differential 
orbit improvement method utilizing observational data is described. This 
method produces a mean set of orbital elements by an iterated least- squares 
fitting of the equations of condition. The results of preliminary applications 
of the orbit generator and differential correction to two artificial satellites of 
the Earth, through use of a high-speed digital electronic computer, is shown 
in tabular and graphical form. 
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ORBITAL PREDICTION AND DIFFERENTIAL CORRECTION 
USING VINTI’ S SPHEROIDAL THEORY FOR ARTIFICIAL SATELLITES 


INTRODUCTION 

The spheroidal method for satellite orbits provides a procedure for calculating the orbit 
of any satellite of an oblate planet, when all forces except those of the primary’s gravitational 
field are neglected. Determining the effect of the oblateness of a planet on the orbit of a satel- 
lite sufficiently near to the planei so that the forces of other bodies may be neglected is one of 
the central problems of satellite astronomy. 

Vinti, in a series of research papers (listed in the references at the conclusion of this 
report), has found a gravitational potential for the exterior of an axially symmetric oblate planet 
which is able to produce an "intermediary reference orbit" accounting for more than 99.5 per- 
cent of the deviation of the Earth’s potential from spherical symmetry. The Vinti potential is a 
very accurate approximation for the Earth’s gravitational potential, which both satisfies Laplace’s 
equation and leads to separability of the Hamilton-Jacobi equation in oblate spheroidal co- 
ordinates, the most appropriate system for an oblate planet. Use of this form for the potential 
reduces the problem of satellite motion to the analytic solution of quartic polynomials and avoids 
the use of perturbation theory entirely in deriving an accurate intermediary orbit. The Vinti 
potential is actually much closer to the empirically accepted one for the Earth than any previ- 
ously used as the starting point of a calculation. In the case of the Earth, the resulting inter- 
mediary orbit reproduces the even zonal harmonics exactly through the second and approxi- 
mately through the fourth. The secular solution can be obtained to arbitrarily high order in the 
second harmonic oblateness parameter, and, by means of rapidly converging infinite series, the 
periodic solution can easily be obtained through second order. The solution holds for all angles 
of inclination (in the case of equatorial or near -equatorial orbits, certain simplifications can be 
made in the equations) and contains no critical inclination or long-periodic terms. For such a 
reference orbit, error can never accumulate because of the exactness of the secular terms. 

This method of solution for unretarded satellite orbits has been adapted for computational 
purposes on a high-speed digital electronic computer primarily by means of the FORTRAN 
programming language. The function of the present paper is to provide the computational pro- 
cedure for determining and correcting an orbit in algorithmic form, adopting algebraic symbols 
consistent with those in Vinti’s papers. A summary of preliminary results utilizing observa- 
tional data from artificial satellites is included. 


INPUT PARAMETERS 

The fundamental physical units employed are those of the canonical Vanguard system. In 
this system, the fundamental unit of length is the Earth’s equatorial radius (taken to be 6378.388 
km), and the fundamental unit of mass is the terrestrial mass (taken to be 5.983 x 10 24 kgm). 
The fundamental unit of time is adjusted so that the Newtonian gravitational constant G is set 
equal to unity; this process yields a value for the Vanguard unit of time of 806.832 seconds. To 
obtain a physical significance for this time, consider a satellite "orbiting" the Earth at its sur- 
face. This time unit is then seen to be the time required for such a satellite to traverse one 
radian. 

The inertial co-ordinate system takes the Earth’s polar axis as the Z-axis (which is also 
the planetary axis of symmetry and the axis of rotation). The X-Y plane is the equatorial plane, 
with the X-axis pointing toward the vernal equinox (the first point of Aries), the Y-axis orthog- 
onally to the east to form a right-handed system, and the Earth’s center of mass at the origin. 
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The following constants are required in the computations: 

GM, where G is the Newtonian gravitational constant and M is the Earth's mass. 
From the preceding remarks, it is seen that - 1 in the Vanguard system. 

c = r e /i7, where r e is the equatorial radius of the Earth (unity in the Vanguard sys- 
tem) and J 2 is the coefficient of the second zonal harmonic in the infinite series 
expansion of the Earth's potential. The value of J 2 is approximately 1.0823 x 10“ 3 . 

t A , the initial time. 

t f , the final time. 

At, the time increment used in generating position and velocity components for equal 
time intervals following t. and preceding t f . 

X., Y., Z., the initial conditions of position. 

X., Y., Z., the initial conditions of velocity. Note that the set X. , Y. , Z. , X. , Y. , Z i of initial 
conditions is also referred to as the set of injection conditions 'if t. 1 , the initial 
time, is taken to be the time of injection of the satellite into orbit. 


CO-ORDINATE CONVERSION 

We now compute the following quantities. The square of the magnitude of the position 
vector: 


r 2 = X 2 + Y 2 + Z 2 . 

1111 

The dot product of the position and velocity vectors: 

r r = X. X. + Y Y. + Z Z . 

ii ii ii ii 

The oblate spheroidal co-ordinates (p, v 9 <P ) and their time derivatives: 

p\ = ^ [(r? “ c 2 ) + /(r 2 - c 2 ) 2 + 4c 2 Z 2 J 


V* ' ~- 2 [- ( r i “* c 2 ) + /(r 2 - c 2 ) 2 + 4 c 2 Z 2 ] . 

Then and t?. are found by extracting square-roots, with the condition that the sign of ^ is 
the same as the sign of z . . 

1 * r i r i ( r i ~ ° 2 ) + Z, Z x 

/9 ‘ 2f> ' ‘ /(r* - c 2 ) 2 + 4c 2 Zj 

1 f . (r* -c 2 ) + 2c 2 Z l Z~ 

V . - r . r + 

2c ^i |_ 1 1 /( r 2 - C 2 ) 2 + 4 c 2 Zj 
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In the above, p i t 0, but if t;. = 0, then v. = z. /p.. 


Y 

sin£. - — - .... 1 . . 

/(Pi + c 2 ) (1 ~ rf) 


cos <p. — , .. 1 .. 

/(P? + C 2 ) (1 - 77?) 

From the above trigonometric relations, we obtain 0 . within the limits 0 < 4>i < 2tt . 


THE JACOBI CONSTANTS OF GENERALIZED MOMENTA 
Compute: 

a, = I (X? + Y 2 + Z»)-M p. (p* + c 2 T7 2 )' 1 
a 3 = Xj Y, - Yj Xj 

a 2 = [(P? + c2, Tf ) 2 V 2 i + a 3 ‘ 2a , C 2 7? 2 (1 ~ r^)J (1 - rp-* /2 

FACTORING THE QUARTICS: PRIME CONSTANTS 
Compute: 



K = c2 Po 2 

A = - 2k oPoYD t _ k o ( 3 x ^Yd ~ 2x i + 4 >] 

B = k o Po C 1 “ Yd> [* + k o Yd ( 4 " ^)] 

b 2 = VB 

3 = Po *£* \} ~ k o x d Yd + k o ^Yd < 3 Yd ~ 2 “ 8 Yd + 4 >] 

p = a' 1 p 2 x^ 2 [l - k 0 y 2 (4 - x 2 ) - 16 kj y 2 (2y 2 - 1) 

- U o4Yd (^Yd " ^ ‘ 2 °Yd + 12 >] 
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S = k 0 Vd ( 3x 2 ■ 4 ) " 16 k o y£ ( 2 Vd " !) 

- k ox£y£ ( 2x ^yD " 5x£ - 28y 2 + 20) 


e = [l - x 2 (1 + g)] 


If (a* - a| ) = 0, then v 0 = 0 and v? = k „ d - k 0 x ^- 


If (a 2 - a 2 ) ^ 0, then calculate v 0 from 




2 "2 ~ 2a l ^ 

2 ( a 2 " a l> 


1 + 


/ ~ 8aj c 2 (a 2 - a 2 ) 

y («2 ■ 2 a i c2 > 2 


Also, 


a 2 -2 a. c 2 
-2_ 2 1 


1 + 


2 ( a 2 - 


8a t c 2 (a 2 -a 2 ) 
(a 2 - 2oj c 2 ) 2 


MUTUAL CONSTANTS 


0 = 


If / 0, then compute: 


Aj = (1 - e 2 ) 172 p 



[( 1 - e 2 ) 172 ] 


where P n (x)is the Legendre polynomial with argument x of degree n, and where R n (x) = x" P n (1/x). 
The infinite series above (and those that follow) is computed by an iterative method, with com- 
putation of terms ceasing when the absolute value of the ratio of successive terms minus unity 
is less than or equal to some pre-selected tolerance, i.e., computation ceases when 


(Vi-. 

(Vi 


< £ 


where e might be 10 _7 . Convergence should be attained by consideration of the first several 
terms, in most cases. To increase computational speed, the first term (for n = 2) of the above 
series may be given explicitly by 


(1 -e 2 ) 172 



If v 0 - 0 (corresponding to an orbit in the equatorial plane), compute instead: 


where the first term (for n = 2) is given explicitly by 


£(l-e 2 )* 72 


& 
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If '; 0 f 0, compute: 


(1 


e 2 ) 1/2 „-l 


l L[i ) p - feh f (i _ e2)1/2 3 


where P n and R are defined as above, and where the first term (for n = 0) of the above series 
is given explicitly by 


(1 - e 2 ) l/2 p 1 . 


If = 0, compute instead: 


A 2 = (1 - eV- p-> 2_' (~) *n [d - e 2 ) 1 ^] 


where the first term (for n = 0) is given explicitly by 

(1 - e 2 ) 1 ' 2 p-». 


Then compute: 

00 

A 3 = d - e») 1/l P' 3 2] »„ R +2 [(1 - e 2 ) ,/2 ] 

n= 0 

where, for r? 0 ^ 0, D n is computed as follows: 



( n an even integer) 


(n an odd integer). 

If v 0 = 0, use instead: 


(n an even integer) 


(n an odd integer). 


D =D„ 




D „= D 2,^ (-!) - m (|) 2 


(4m)! 

[( 2 m )!] 2 


° n °2,, = r (-i 

4^ W [(2m ♦ 1)!] • 
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Then compute: 


B. = I + A q 2 + JJL q 4 
1 2 16 128 


b 2 = i q 4 


4 64 


b 3 = i - (l -v- 2 2 V 1/2 - y m v 2 2m 

m = 2 

where 


m- 1 

(2m)! yi 

2 2m (m!) 2 4f 


(2n)! 
2 2n (n! ) 2 


V 


2n 

0 


A n = | (1 - e 2 ) 1/2 p ‘ 3 e (- 2 bjb 2 p + b^) 

A j 2 = C 1 ” e 2 ) 1/2 p -3 b 4 e 2 

If (b,/b 2 ) < 1, then compute: 

A 2l = C 1 " e 2 ) 1/2 p _1 e jbj p _1 + (3 b 2 - b 2 ) p ' 2 

-| b l b 2 t 1 + 5 e2 ) P' 3 + f b 2p-4 (4 + 3e 2 )J 
A 22 = (1 - e 2 ) 1/2 p ' 1 |j- e 2 (3 b 2 - b 2 ) p * 2 

- | e 2 bjb 2 p ' 3 + b 4 P ' 4 ( 6e2 + e4 )J 

A 23 = i (1 - e 2 ) 1/2 p _1 e 3 (- bjb 2 p ' 3 + b 4 p* 4 ) 

O 

^ = tL (1 ' e2)1/2 P " 5 b 2 p4 
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A - (1 - e 2 ) 1 72 p' 3 e 


2 + b j P _1 ^3 + j e^j - p~ 2 b 2 + c 2 ^ (4 + 3e 2 ) 


A„ = (1 ~ e 2 ) l/2 p 3 


— e 2 + 

4 


fbiP" 1 e 2 - p" 2 ^-b 2 + c2 )(|e 2 +|e 4 j 


A 33 “ ( 1 “ e 2 ) l/2 P~ 3 e 3 


\u pl b i-l p ' 2 (i b2 + c2 )_ 


A 34 ’ “ 32 (1 “ e2)1 2 e4 pS 


(f b2 + c2 ) 


Ii ( !) i /i\ ) -- 1 (corresponding to a near-equatorial orbit or an equatorial orbit), then compute 
instead: 

A 2J " - c 2 ) 1 2 p" 1 e Jbj p' 1 + (3 - h\ b" 2 ) c 4 p -4 (1 - t 7 2 ) 2 J 

A 22 = J (*■ ~ e 2 ) 12 P _1 e2 (3 " b 2 bj 2 ) c 4 p“ 4 (1 - t] 2 ) 2 


A 23 and A 24 as given above. 


A 31 - (1 - e 2 ) 172 p‘ 3 e 




' 2 p~ 2 (1 - V 2 ) - (4 + 3 e 2 ) c 2 p 


A 32 ” (1 “ e 2 ) 1 2 p " 3 e 


2x1 2 ^-3 „2 


I + | c 2 p- 2 (1 - T, 2 ) - e 2 + c 2 p- 2 ] 


A 33 =1(1 + e 2 ) 1 2 p-» e 3 c 2 




-■] 


A 34 — — d-eV V 5 


32 


Now compute: 


If (bj/b ) < 1, compute: 


-1 


277! "l = (- 2a, ) ] 2 (a + b, + A, + c 2 rj* B, B/ 1 ) 


2vv 2 = (a 2 - a 2 ) 1 2 tj‘ ! A 2 B" 1 (a + b, + A, + c 2 t? 2 A 2 Bj B" 1 )' 
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If (b x /b 2 ) > 1, compute instead: 

^ 1JV 2 - a 2 - 2 a 1 a 2 2c2 ( 1 “ vl)l 1/2 A ^' 1 ( a + b ! + A t + c 2 ^ 2 A 2 B' 1 )* 
Then compute: 

e' = ae(a + b 1 )" 1 . 

Note that the condition e' < e must be fulfilled. 

THE JACOBI CONSTANTS OF GENERALIZED CO-ORDINATES 
If e ^ 0, then compute: 

Pi (p\ + c 2 ) 

sin E. = ■/: . ■- ... . — ■ - 

1 ae /(-2a 1 )(p2 ~ A/D . + B) 

cos E. = (1 - p i a -1 ) e" 1 . 

From the above, we obtain E { within the limits 0 < E i < 2tt. 

(1 - e 2 ) 1/2 sin E l 

sin v. ~ - 

1 - e cos E A 


cos v i 


cos E i - e 
1 - e cos 


From the above, we obtain v. within the limits 0 <_v. < 2v. 


If e = 0, then: 


and 


If 77 / 0, then compute: 



COS Ipy 


Vj, (P? + rf C 2 ) 

/(- 2 a i)( 7 )? - vf) 


sin = 


From the above, we obtain <p. within the limits 0 < i/<. <2v. 
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(1 - T?g) 1/2 sin 0. 
1 A ~ ^ sin 2 i/T 


V: = 


cos x. 


cos \p. 


' A - r£ Sin 2 i//. 

From the above, we obtain x , within the limits 0 < x . < 2 tt . 

If v 0 = °» then: ^ = 4>. and x . = <*>. . 

Now compute: 

sin nvj for n = 2, 3, 4- 
sin m//. for n = 2, 4. 

If (bj/bj) < 1, then compute: 

jB t = (- 2 a,)~ 1/2 [bjEj + a(Ej - e sin Ej) + A, v. 

+ A,, sin v. + Ajj sin 2v.] + c 2 (a 2 - a 2 )~ 1/2 rj* |b, 0, 

- — (2 + q 2 ) sin 2< />, + — q 2 sin 40. - t. 

8 64 J 

/S 2 = - a 2 (-2a,)-" 2 [Aj Vj + A 2I sinv. + A J2 sin 2v. 

+ A 23 sin 3v. + A 24 sin 4v,] + (a 2 - a 2 )* 1/2 V 0 a 2 |b 2 

- J_ q 2 (4 + 3q 2 ) sin 2 0, + — — q 4 sin 40, • 

32 256 J 

If (b, /b 2 ) > 1, then compute instead: 

i 8, = (- 2a,)’ 1/2 [b, E, + a(E, - e sin E,) + A, v. 

+ Ajj sinv. + A, 2 sin 2v.] + c 2 tj 2 a* 1 [1 - 2a, a' 2 c 2 

- i (2 + q 2 ) sin 20j + — q 2 sin 40, - t. 

8 64 J 


(t ~ 
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/3 2 = - a 2 ( - 2a !)“ 1/2 C A 2 v i + Aj, sin Vj + A 22 sin 2v. 

+ A 23 sin 3vj + A 24 sin 4v. ] + [1 - 2a 1 a 2 2 c 2 (1 - i ) 2 )]' ,/2 

~J2 q2 (4 + 3q2) Sln + 2^6 ^ Sin 4 ^] ' 


If 7j ^ 0 and if (b t /b 2 ) < 1, then compute: 


p 3 = ^ - a 3 ( a 2 - a 3 y in 


+ B 3^i +12 ^ 2 ^i 


(1 -i?g)- 1/2 (i - V~ 2 2 )~ 1/2 X. 

+ c 2 a 3 (- 2aj)* 1/2 [A 3 Vj 


+ A 31 sin V. + A 32 sin 2v. + A 33 sin 3v. + A 34 sin 4v.] . 


If rj £ 0 and if (b t /b 2 ) > 1, compute instead: 

- <P L - a 3 a 2 1 Cl *“ 2a i a 2 2 c2 ^l ~ ^o^" 172 C(i “ Vl )' l/2 (1 “ ^V 172 ^ 

+ B 3 0.] + c 2 a 3 (~2a ] Y l/2 [A 3 v. + A 31 s i n + A 32 s in 2v. 


+ A 33 sin 3v. + A 34 sin 4v.] 

If v = 0, compute instead: 

Aj = - $2 ( s e n a 3 ) “ a 3 (‘ 2a i)" 172 CAj V. 

+ A 21 sin v. + A 22 sin 2v.] + c 2 a 3 (- 2a 1 )" 1/2 [A 3 v. 

+ A 31 sin v. + A 32 sin 2v. + A 33 sin 3v. + A 34 sin 4v.] 


, 3 

where sgn a - - — - . 

Kl 

THE ORBIT GENERATOR OF POSITION AND VELOCITY COMPONENTS 

In this section, parameters arise which are time-dependent. Initially, the value for time 

t is equal to t., but on subsequent iterations t = t. + n (At), n = 1, 2, 3, Here At is the 

time increment input parameter used in generating position and velocity components for equal 
time intervals following t and preceding or coincident with the final time t f . 

If v 0 t 0 and if (b,/b 2 ) < 1, then compute: 


M s = 27 TVj (t + /3j - c 2 AjO." 1 r? 2 Bj B" 1 ) 

4> s = 2 ttv 2 [t + + ^aj 1 A' 1 (a + b t + Aj)] 
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If t) 0 t 0 and if (b x /b 2 ) > 1, then compute instead: 

B 2 (* + ^l) “ c2 ^ 2 a ] i 1 Vi Bj “ 

(a + b x + Aj)B 2 + c 2 tj 2 AjB^ 

^ as given above. 

If t) 0 - 0, then compute instead: 

M s = (- 2* 1 ) 1/2 (t + /3j)(a + b x + Aj)* 1 

= (1 - 2a 1 a' 2 c 2 ) 1/2 [/3 2 + a^ft + / 3 X )( a + bj + Aj)” 1 ] 

We now solve the following equation for (m s + E 0 ): 

M s + E 0 - e' sin (M s + E Q ) = M s . 

If we let 8 = M s + , then we can solve this equation (known as Kepler's equation) by use of the 

iterative Newton-Raphson method. 

g _ g (?„- e* sin 8 n -M s ) 
n+1 " (l-e'cos8 n ) 


M s = (- 2a ,)*'* 


(8 n - e' sin 8 n -M s ) 2 (e'sin£ o ) 
2(1 - e' cos 8 n ) 3 


For the initial value, (8 ) 

7 v n y n=0 


M s . Iteration ceases when 



< 


e 


where e is a pre-selected tolerance (e.g.,10~ 7 ). Convergence should be attained with several iterations. 
Now use the anomaly connections: 

cos v' = (cos 8 - e)(l - e cos 8)" 1 

sin v' = (1 - e 2 ) 12 (1 - e cos 8)" 1 sin 8. 

From the above, we obtain v* within the limits 0 < v' < 2n . The angle v' is then placed within 

the same circle of revolution as the angle 8 (which is not taken modulo 2 t7). 

Then compute: v Q = v' - M s 
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If v 0 £ 0 and if ( b i A 2 ) < 1, then compute: 

P 0 = (-2 a ,)' 1/2 ( a 2 - a *) 1 ' 2 tj " 1 A, B’ 1 v Q 

M t = (a + b,) _l [-(A, + c 2 tj 2 A 2 B 1 B" 1 )v„ 

+ I C 2 (-2a,) 1/2 (a 2 - a 2 ) ' 1/2 t)\ sin (2 yp t + 2 -/-„)] 

If v 0 / 0 and if (b,/b 2 ) > 1, then compute instead: 

r -i 1/2 

P o=(- 2a i)" l/2 [l - 2a,a 2 2 c 2 (1 - t] 2 )J cij AjB’ 1 v q 
M, = (a tb,)* 1 [- (A, + c 2 7j 2 A 2 Bj BJ *) v 0 

+ Ic 2 (-2a,) ,/2 [l - 2a, a' 2 c 2 (1 - i?|)J ' 1/2 v Jaj 1 sin (2</» s + 2 </-„)] 

If 7j 0 = 0, then compute instead: 

= (-2 a,)' 1/2 (1 - 2a, a' 2 c 2 ) 1/2 a 2 A 2 v 0 
M, = - (a + b,)' 1 A, v fl 

Continuing, if v 0 / 0 and (b,/b 2 ) < 1, or if y 0 = 0, then compute: 

E, = (1 - e' cosg)" 1 M,-ie' (1-e' cos8) _3 M 2 sin 6 

If v 0 i 0 and if (b,/b 2 ) > 1, then compute instead: 

E, = (1 - e' cos 8)“ l M, 

Now use the anomaly connections again: 

cos v" = [cos (g + E,) - e] [1-e cos (£ + E,)] ” 1 

sin v" = (1 - e 2 ) 1/2 [1-e cos (g + E,) ] ' 1 sin (g + E,) 

From the above, find the angle v " and place it within the same circle of revolution as the 
angle (8 + E, ). 

Then compute: v, = v" - v' 

li v 0 ^ 0 and if (b,/b 2 ) <1, then compute: 

0, = (-2a,)- 3/2 (a 2 - a 2 ) 1/2 rj~ 0 l Bj 1 (A 2 v, + A 2 , sin v' 

+ A 22 sin 2v' ) + i q 2 B' 1 sin (2^ s + 2<p 0 ) 
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If t? 0 t 0 and if (bj/b 2 ) > 1, then compute instead: 

= (-2a J )‘ ,/2 [1 - 2<Xj a* 2 c 2 (1 - 7 j 2 )] ,/2 a 2 Bj l (A 2 v, 

+ Ajj sin v' + Aj 2 sin 2 v' ) + A q 2 B~ 1 sin (2 >/>, +2 1 p Q ) 

If tj 0 = 0, then compute instead: 

i/ij = (-2aj)" 1/2 (1 -2», a* 2 c 2 ) 1 -' 2 a 2 (A 2 Vj + A,j sin v' + Aj 2 sin 2v' ) 

Now if (bj/b 2 ) < 1, we continue this procedure one step further to obtain terms of second order, 
as follows. Compute: 

M 2 = - (a + bj) -1 {Aj Vj + A„ sin v' + A, 2 sin 2v' 

+ c 2 (- 2 a.j) 1/2 (a 2 - a 2 )‘ 1/2 V § Jb, -A'/', cos (2</«, + 2*p 0 ) 

-Aq 2 sin (2i/< s + 2 </>„) q 2 sin (4 + 4^ 0 )j| 

E 2 = [1 -e' cos <8 + Ej ) ] ~ 1 M 2 


Let 


E = 8 + Ej + E 2 


and use the anomaly connections once again: 

cos v'" = (cos E-e) (1 - e cos E)’ 1 
sin v'" = (1 - e 2 ) 172 (1 ~ e cos E)' 1 sin E 

From the above, find the angle v" and place it within the same circle of revolution as the angle E. 
Then compute: v 2 = v w - v" 

+2 = (- 2 “.)' ,/2 (“ 2 2 - “^) 1/2 V B-1 (A 2 v 2 + A 21 v, cos V ' 

+ 2 A 22 Vj cos 2 v* + A 23 sin 3v' + A 24 sin 4 v' ) 

+ Iq 2 B* 1 ^cos (2^ s +2 4 > 0 ) + -| q 2 sin (2 </>, + 2</< 0 ) - q 2 sin (40 s + 4 </.' 0 )J 


Finally, let: 


= M s + v 0 + v, + v 2 


'P ~ P t + 'Po + ^1 + Pi 
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Now if (b t /b 2 ) > 1, we omit computation of M 2 , E 2 ,v 2 , and </< 2 . In such case, these terms be- 
come of the third order and hence negligible. 

Instead, we let: 

E = £ + Ej 
v =M, + v 0 + Vj 

s + -Ao + '/'i , 

Continuing, compute: 

sin x = (1 - Vq) 1/2 (1 - V q sin 2 1/2 sin 

cos y = (1 - rj 2 sin 2 0)“ 1/2 cos i p 

From the above, find the angle y and place it within the same circle of revolution as the 
angle 4> . 

If (bj/bj) < 1, then compute: 

p - (1 + e cos v)“ 1 p 

If (bj/b 2 )> 1, then compute instead: 

p - a (1 - e cos E) 

Then, if v 0 / 0 and if (bj /b 2 ) < l, compute: 

V°V 0 sin 4> 

4> = 0 3 + a 3 (a 2 - a 2 )' 1/2 v 0 |(i - vl)' U2 (1 * V~ 2 2 )~ 1/2 V 
? 2 t; 2 4 sin 2^j - c 2 a 3 (-2 0 ^)' 1/2 (A 3 v 

+ A 31 sin v + A 32 sin 2 v + A 33 sin 3 v + A 34 sin 4 v) 

If ^ 0 and if > 1, compute instead: 

v as given above. 

0 = /3 3 + a 3 a* 1 [l -2a lQ - 2 c 2 (1 - r, 2 )] ' 1/2 [(1 - ^|)“ 1/2 (1 - ^ 2 )“ 1/2 y 

+ B 3 - c 2 a 3 (- 2 oij)" 1/2 (A 3 v + A 3 j sin v 
+ A 32 sin 2v + A 33 sin 3v + A 34 sin 4v) 
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If v Q = 0, compute instead 


v = o 

^ > = ' fi 3 +^2 ( s B na 3 ) + a 3 (-2a 1 )-»/2 (AjV 

+ A 21 sin V + A 22 sin 2 v) - c 2 a 3 (-2a 1 )‘ 1/J (A 3 v 

+ A 31 sin v + A 32 sin 2v + A 33 sin 3v + A 34 sin 4v) 

a 

where sgn a = 

I a 3 I * 

The oblate spheroidal co-ordinates must satisfy the following conditions: 

p > 0 

- 1 < v < + 1 

Now the co-ordinates and velocities may be found as follows: 

h l = ( p 2 + c2 ) (1 - v 2 ) 

p = a e (- 2 a 2 ) 1/2 (p 2 + Ap + B) 1/2 (p 2 + t? 2 c 2 )" 1 sin E 
V = cj] 0 (- 2a 1 ) 1/2 (i 7 2 - 7? 2 ) 1/2 (p 2 + t] 2 c 2 ) -1 cos 0 
X = /(p 2 + c 2 ) (1 - 7 } 2 ) cos 0 

Y- /(p 2 + c 2 ) (1 - t? 2 ) sin 

Z = p V 



Z = p T) + V p 


This completes the algorithm for predicting orthogonal position and velocity components of the 
satellite based upon a set of initial conditions. 
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COMPUTATION OF DIRECTION COSINES 


Often the set of initial conditions (position and velocity components) provided is only ap- 
proximate at best, and thus the orbit that is predicted based on these initial conditions will 
similarly contain inaccuracies. In order to remove these inaccuracies and to account for the 
effects of forces not considered by the analytical theory (e.g., aerodynamic drag, solar radia- 
tion, meteoric bombardment, etc.) the orbital parameters are corrected by comparison with 
those found by direct observation. The orbit improvement method produces a mean set of 
orbital elements through an iterated least- squares fitting of the differential solution to numer- 
ous observational values. 

To perform the differential correction process, the following data must be available in 
addition to the constants listed in the section on ’’Input Parameters” above: 

f a l _ r /r e , the flattening coefficient of the Earth (where r p is the polar radius of the 
Earth, and r e P is the equatorial radius), taken to be approximately 1/298.3 = 3.3523299 x 10 ” 3 . 

oj , the rotational rate of the Earth in radians per mean solar hour (taken to be 0.26251614). 

(x^). , i = l,2 ,...,s, the geodetic (or geographic) longitude of the terrestrial tracking 
stations in radians, as measured eastward from Greenwich (a negative sign must be prefixed 
if measured westward from Greenwich). We assume that there are s tracking stations report- 
ing observational data used for comparison purposes. 

(( 9 )., i = 1, 2, . . . , s, the geodetic (or geographic) latitude of the stations in radians, 
measured as positive north of the equator and as negative south of the equator (-tt -/ 2 +77 / 2 )* 

(h) . , i = 1, 2, . . . , s, the altitude of the stations in feet, measured as positive above sea 
level and as negative below sea level. 

(\ ) d , d= 1, 2, . . . , the angle in radians, measured eastward from the vernal equinox (the 
first point of Aries) to the Greenwich meridian at midnight Greenwich mean time for each day d 
during the period that observations are provided. The apparent sidereal time (the hour angle of 
the first point of Aries) at midnight Greenwich mean time for each day throughout the year is 
tabulated in ’’The American Ephemeris and Nautical Almanac.” 

t , a reference time preceding or coinciding with the time of the first observation pro- 
vided, which is used as the zero point in determining t , the relative observation time. It may 
be the time of injection if the tracking data includes observations made during the first several 
orbits of the satellite. The purpose of determining a relative observation time t is to elimi- 
nate any reference to the calendar. 

We now describe the observation data cards, which are effectively input for the differ- 
ential correction scheme. There are several methods of recording satellite tracking data; 
we present here one of the most common methods, referred to as Minitrack observation data. 
(Refer to the appendix of this report for discussion of another method.) A set of observation 
data of this type includes the following parameters for each recorded spacecraft observation: 

t' , the date and time of observation. As given, t' is a calendar time. We remove any 
dependence on the calendar by determining t = t' - t 0 , where t is the relative observation 
time and t Q is a reference calendar time. Then t becomes a time interval, measured in Van- 
guard units of time from the zero point t Q . As mentioned above, t Q is chosen so that for all 
observations t > 0. It is convenient to have the choice of coincide with that corresponding 
to the initial position and velocity conditions X i? Y., Z., X., Y., Z. . When this choice is made, 
then t 0 is known as an initial or epoch time. 
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k , the code number for the tracking station reporting the observation. Generally, the 
range of k is the set of integers 1, 2, 3, . . . , s . 

L 0 , the observed direction cosine in the X-direction. 

M 0 , the observed direction cosine in the Y-direction. 

w andw^, the weighting factors corresponding to observations L 0 and M 0 f respectively. 
This information is optional; if not provided, then it is assumed that w L and are each unity. 

The co-ordinate system employed for the observation data is centered at the tracking 
station on the Earth’s surface, with the X-Y plane tangent to the surface. It is a right-handed, 
orthogonal system with the X-axis extending in an easterly direction along the line of latitude, 
the Y-axis extending in a northerly direction along the line of longitude, and the Z-axis normal 
to the surface and pointing toward the geodetic zenith. 

We first compute the so-called "auxiliary functions" S and C (refer to the "Explanatory 
Supplement to the Astronomical Ephemeris and the American Ephemeris and Nautical Alma- 
nac") from the relations: 


C = [cos 2 & B + (1 - f) 2 sin 2 


S = (1 - f) 2 C 


Here and in the following, we eliminate use of the subscript "i" referring to an individual one 
of the s tracking stations, and assume that the computations given are performed for each re- 
spective station. The value of H is then converted from its input units of feet to units of the 
Earth’s equatorial radius (the conversion factor is 4.77865 x 10‘ 8 ) so as to conform to the 
canonical Vanguard system of units used throughout (see above under "Input Parameters"). 
Then the geocentric latitude is given by: 



Now the geocentric distance of the observation point (i.e., tracking station), in units of the 
Earth's equatorial radius, is found: 

p - [(S + H) 2 sin 2 $ D + (C + H) 2 cos 2 1/2 - 

The angle S , between the vernal equinox and the observation meridian plane, is computed in 
radian measure by the following expression: 


5 - <*») d +"<AT) + V 


Here, (\ 0 ) d is as defined above with the value chosen (designated by the subscript "d") 
corresponding to the midnight immediately preceding observation time. Also, at is the com- 
puted time, in hours, between observation time and midnight preceding observation time. Thus, 
the second term in the expression for 8 accounts for the fact that the Greenwich meridian ro- 
tates while the vernal equinox remains fixed in inertial space. 

The inertial geocentric co-ordinates of the observation point are now converted from a 
spherical to a Cartesian system by means of the following equations: 
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X T = p cos 8 C cos 8 
Y t = p cos 8 q s in 8 
Z T = p sin 8 q 


The angle 4 > x , measured in the observation latitude plane between the vernal equinox and the 
tracking station’s X-co-ordinate axis, is given in radian measure by: 


<p r 8 + — 

2 


The local co-ordinates of the satellite (X M ,Y M , Z M ) can now be determined from the iner- 
tial position vector (X,Y ,Z ) computed by the orbit generator. Here "local" refers to co- 
ordinates measured at the tracking station. The orbit generator will produce the position com- 
ponents (x, Y,z) at the observation time. Then the local co-ordinates are given by the matrix 
relation: 


Y„ = 0 sin 6> cos 


0 - cos 


cost// sini// 0 


sin i p cost// 0 


0 1 


Here, the difference of column matrices on the extreme right represents a translation from 
the Earth's center to the tracking station position; the center matrix on the right represents 
a rotation in the latitude plane about the polar axis through an angle of 'P x to bring the inertial 
X-axis into coincidence with the station’s x M -axis; the left matrix represents a rotation in the 
longitude plane about the X M -axis through an angle of (rr/2 - ) to bring the inertial Z-axis into 

coincidence with the station's Z M -axis. This matrix equation is equivalent to: 


-sin i p % sin 8 D cos ip x sin cos Y - Y T 

sin i/i x cos -cos cos ^ sin Z-Z T 


or, explicitly stated: 

X M = (X - X T ) cos p x + (Y - Y t ) sin ^ 

Y m =-(X-X T ')sin ^ sin 8 0 + (Y - Y ? ) cos <P x sin + (Z - Z ? ) cos & D 
Z„ □ (X - X T ) sin <A x cos 6 d - (Y - Y ? ) cos P x cos 6^ + (Z - Z T ) sin & D 

We now find the computed values of the direction cosines L c (in the X-direction) and M c ( i n the 
Y-direction) in terms of the local co-ordinates. 

X M 

Lc ’ c*2 " *1 + ZS) 175 


(X’+Y’+Z’)"’ 
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Of course, the computed value of the third direction cosine N c (in the Z-direction) is pre- 
determined by L c and M c through the relation: 

N c =(1 -L 2 -M 2 )>' 2 


THE STANDARD DEVIATION OF FIT 

The differences between the observed and computed values of the direction cosines can 
now be found: 

AL = L 0 - L c 
AH=M 0 -M C 

These differences are sometimes referred to as "residuals", although this term is also used 
in a different sense in the method of fitting by least squares. We compute these differences 
for each observation in the set of observation data. The number of observations in the set is 
variable, and it may be determined by an input parameter n . 

The average residual is given by: 



n 



i * 1 


+ am.), 


where the subscript "i” ranges over individual observations. 

The standard deviation of the residuals from their mean value is found from: 


f n 

£ t AL i + < AM i -® )2 ] 


i = 1 


The standard deviation of the residuals (from zero) is called the standard deviation of fit, 
and is given by: 


cr 


f 



[(AL ^ 2 + (AMj) 2 ] 


As is customary, the larger multiplicative factor (2n - 6)~* is used to indicate the excess 
of simultaneous equations of condition over the number of independent coefficients (see below 
under section titled, "Fitting by Method of Least Squares"). 

We may also determine an acceptable range of values for the residuals, bounded by an 
lower limit and an upper limit r 2 , based upon the standard deviation. If a residual falls 
outside this range, it may be rejected, with statistical validity, from inclusion in the fitting 
process. For example, we may choose: 


r 1 ^R-jo- 

r 2 = R + j cr 
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For normal (Gaussian) distributions, 68.27 percent of the cases are included within one stand- 
ard deviation on either side of the mean ( j = 1 above), 95.45 percent of the cases are included 
within two standard deviations ( j = 2 above), and 99.73 percent of the cases are included within 
three standard deviations ( j = 3 above). For moderately skewed distributions, the above per- 
centages may hold approximately. If certain of the residuals are rejected on this statistical 
basis, the standard deviation of the accepted residuals only may be computed as a "working” 
standard deviation of fit. Its value is computed exactly as is above, with certain terms 
omitted in the summation, and should be substantially smaller in magnitude than °- f . 


ANALYTICAL PROCEDURE OF DIFFERENTIAL CORRECTION 

The first-order Taylor series expansion of the equations of condition may be written: 



i=j 


where q ; ( i = 1, 2, . . . , 6) are the mean or Izsak elements given below. 

Qj = a, the semi- major axis. 
q 2 = e, the eccentricity. 

q 3 = v 0 3 sin I, where I is the inclination of the orbital plane to the equator. 

q 4 = Z 3 ! , corresponds to the negative of the time of passage through perigee in Keplerian 
motion. 

q s = P 2 ! corresponds to the argument of perigee in Keplerian motion. 

q 6 = p 3 , corresponds to the right ascension of the ascending node in Keplerian motion. 

We may expand the above partial derivatives by the chain rule as follows: 

3q, BX m 3q. 3 Y m 3q. 3Z M 3q. 

Bq s " BX m Bq. 3Y„ Bq ; 3Z„ Bq. 

From the equations for L c and M c in terms of the local co-ordinates given earlier (refer to the 
section titled, "Computation of Direction Cosines"), we find directly: 

^ = <*2 + y m +%r 1/2 - K (K + y m +z 2 m r 3/2 
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! r =- v . < x « + y m 2 + z «)" 3/2 

|-V.(X^V^Z2)' 3/J 

BUI BL 

Wx$ + y 2 + z 2 )- 3/2 = ^ 

BM 

= (Xi + Y 2 + Z 2 U )- 1/2 - V 2 (X 2 + Y 2 + Z 2 u )- 3/2 


BM 

c 


-Y m Z m (X 2 + Y 2 + z*r 


i /2 


Since the co-ordinates X t ,y t , Z T and the angles V>„ and 0 D are independent of orbital param- 
eters (and merely geodesic functions), we have the matrix relation: 


5*" 


cos </* 

sin 

0 


BX 

BY„ 

Bq s 

= 

-sin0 x sin6? D 

cos sin (9p 

cos # D 


BY 



s i n cos 

“ COS 0 X COS t? D 

sin d D 


BZ ' 

Bq, 


t . , Bx BY BZ 
We find — , _ , — 
Bq ; B q Bq 


by substituting: 


p - a(l - e cos K) and rj □ 77 s i n y 


in the relations: 


X = /(p 2 + c 2 ) (1 - 7 j 2 ) cos 0 
Y = )^p 2 + c 2 ) O - T} 2 ) sin# 

2 - prq 


and then determining 3E/aq. , ty/d q. , and3#/3q.. Here E and ^ are unif or mizing variables, 
analogous respectively to the eccentric anomaly and the argument of latitude in elliptic motion. 
The parameter # is the third oblate spheroidal co-ordinate, the geocentric right ascension. 
The procedure for determining the eighteen partial derivatives BE/Bq, , d^/dq. , and ^ / q, is 
a rather lengthy one which is initiated in the next section. 
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Before embarking upon this procedure, the following comments are in order. The ana- 
lytical partial derivatives in the differential correction given in this report correspond to the 
case where (bj/b 2 ) < 1 only. This excludes equatorial and near- equatorial orbits. More 
specifically, orbits where the inclination I is such that 

0 < I <I C or 180° - I c < I_< 180°, 

where i c might be as large as 1° 54* for an orbit sufficiently close to the Earth, are excluded. 
The analytical partial derivatives for the case in which (bj/b 2 ) £ 1 are, in fact, simpler in 
form, and they are derived in an analogous manner from the equations for this special case 
presented earlier in this report. 

The differential correction process may be carried through to terms of second order, or 
it may be simplified to omit terms of purely second order. This is not quite the same as carry- 
ing the process through to terms of first order, since some second-order effects are included 
even when terms of purely second order are dropped. Generally, the speed of computation in 
the differential correction will be increased considerably by neglecting purely second -order 
terms, without risking any great loss in precision of the final differential coefficients for the 
conditional equations. It should be remembered, however, that even a slight loss in the pre- 
cision of the differential coefficients may cause an additional iteration in the least-squares 
fitting to be necessary. An option to choose the method desired in the differential correction 
may be provided by inclusion of an input variable assuming either one of two values as appro- 
priate for the choice. The terms that are to be omitted in the simplified version will be indi- 
cated as such in the sections that follow. 


PRIME CONSTANTS II 

The following parameters are utilized extensively throughout the differential correction 
process and must therefore be evaluated beforehand. In many cases, the parameters are those 
that were computed previously by approximation methods, and are here re- determined by more 
accurate expressions. 


p = a (1 - e 2 ) 


D = (ap - c 2 ) (ap - c 2 tj 2 ) + 4 a 2 c 2 tj 2 
D' = D + 4 a 2 c 2 (1 - t) 2 ) 

A= - 2 ac 2 D' 1 (ap - c 2 ?7 2 ) (1 - tj 2 ) 

B = c 2 7j 2 ET 1 D' 



b 2 = Vb 

a X = " J M (a + b,) -1 

a 2 = [/x (a + bj) _1 ] 1/2 [ap D' D _1 - c 2 (1 - 17^)]* 2 
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a 3 = a 2 {* ' c2 -nl UpD’ D* 1 - c 2 (l - 7^)] _1 | 1/2 (1 - r^) 1/2 
T}~ 2 - c 2 D (ap D') -1 

q = W2 1 
k = c 2 p -2 

DIFFERENTIAL CORRECTION: TIME -INDEPENDENT PARTIAL DERIVATIVES 
Compute the following in the order indicated: 

= 1 - e 2 
Ba 

l£=-2ae 

Be 


= - (1 - e 2 ) p“ 2 



= 2 a -1 e (1 - e 2 )~ 2 


BD 

Ba 


8ac 2 T7 2 + 2p [2ap - c 2 (1 + t? 2 )] 


BD 

Be 


[2ap - c 2 (1 + 77 2 )] a |E 
oe 


BP 


~ 2 (ap-c 2 ) c 2 j) 0 


BD* 

Ba 


BD 

Ba 


+ 8ac 2 (1 - t? 2 ) 


BD' _ 3D 

Be Be 


bd; 

^Vn 


cD 

dr k, 


- 8 a 2 c 2 T) 
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3D- 1 _ _ d - 2 3D 
Ba Ba 


= -D" 2 

Be 


BP" 1 p , 2 3D 

3t >o B7 ?0 


3b, r , -i 

7“ =b, [— 

Ba L a P " T 7 ‘ Ba 


Bb 


i_ 

Be 


ac 2 (1 - 7) 2 ) 


aD_l + (ap - c 2 7j 2 ) 


3D- 1 

3e 


] 


2 r 

^' ac l 


(ap- c 2 r, 2 ) - 2c 2 r, 0 D- 1 ] (1 - r, 2 ) - 2 r, 0 D" 1 (ap-c 2 7 , 2 )| 


3b, 


2 

Ba 


Bb 2 


2E! + d- 

3a 

BD- 1 

3a 

i51 + D ' 

Be 

3D-‘\ 

3e J 

™ + D' 

3D-i^ 

3t)o/ 


Ba i 1 , 

— — = — u. (a + b, 

3a 2 ^ v 


i } ' 2 { l + bt) 


1 1 ^ b , 

1 _ i . . l \_2 1 


Ba 

"B^ 2 


_L = _ M(a + bl) -_^. 


Ba, 1 _ 

‘ = ^ (a + b t )- 2 

H 2 


3^ 

3t ?o 


V5‘‘ l '* b 1 ) '' 1 “ * a -( D '‘ lr tD 'irr)_ 

-(a +b 1 )- 1 [apD'D- 1 - c 2 (1 - 7 ^)]i/a ^1 + 


[apD'D* 1 - c 2 (1 - t, 2 )] 


2M- 1/2 
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-P- =1 [ M (a +b,)- 1 ] 1/2 JaDT)- 1 |E+ap 
oe 2 1 Be 


| aD*D-' |E + ap ^D" 1 + D' ^p D'lT 1 - c 2 (1 - t; 2 ) 


- (a + b,)-» [apD'D- 1 - c 2 (1 - 17 2 )] 1/2 -I 


K = 5 6* <- + 


D" 1 + D' + 2 c 2 7,^j [ap DT)- 1 - c 2 (1 - t, 2 ) 


- (a + bj)- 1 [apD'D* 1 - c 2 ( 1 - t? 2 )] 


Bb. 

1/2 1 


-2 3a 2 


-2 js 

°~ 2 Be 


3a 0 | Ba apD'D- 1 - c 2 (1 - P) 


c2 ^ a 2 j c _^i 

2 [apD'D -1 -c 2 ( 1 -t^)] 2 _ apD'D -1 -e 2 (l -tj 2 ) 


. f , BD' BD~ 1 

2 pD' D* 1 4 ap D' 1 -e- 4D' /— 
\ a a Ba 


Ba 2 c2t ?o 

M - t) 2 )1 /2 <— 1“ ° 

V V |3e apD'D->-c 2 (l-^) 


2[apD'D- 1 -c 2 (l-7j 2 )P |_ apD'D-‘-c J (l -tj 2 ) 


a D D 1 + ap D — — 4 D 

Be V Be <’<' 
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~| 1/2 


3 a. 




3 ct 2 


i - 


c 2 V 2 0 


ap D'D * 1 - c 2 (1 - 77 2 ) 


C 2 77 a 

( '0 2 

* l 

o 

to 

cFL 
1 

-1/2 

2 [apD'D* 1 -c^I-t; 2 )] 2 

apD'D* 1 -c 2 (l -ijg) 



2 (c 2 - ap D' D *) + ap 77 


A i 3D # , 3D 

D 1 — + D' — 

V dT >o a7 )< 


-1 

'0 


- ^0 V 1 - ^ o )' 1 

^2 1 -2 

■aT = T pc *2 


aD' ^ 


D ‘ 1 l 2D ' + a irj +aD ' 3a 


3 77 , 


f 2 1 


.-2 -1 


= o ac ^ 




D“ J D' ^L +P 

de 


S)*0^' 


-f " ~-a p c ~ 2 rjJ* (d->|21 + D' 3D_1 


<H 2 


a? 7o 3r ?o 


3 q .2 dr >2 

37 = " 


3 q 

3e 


- ^ 




Be* 

3 a 


(a + bj )~ 2 e 


3b, \ 
3 3a ) 


de 1 
3e 


= (a + bj )“ 2 


a 




de ' 

?i7 ?n 


- ~ (a + b,) 


-2 




3 Bj 

~ i 7 



3 q 

3a 



3 q 

3 e 
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3B, 

-(A q 

+ il n 3'l 

Bq 

9t ?o 

U Q 

32 q J 


BB2 

(i q + 

-9-^ 

3q_ 

3a 

u 

16 Q J 

Ba 

3B 2 

~bT = 

G" 


Bq 

Be 

3B 2 

_/J_ q 

+ 9 2 

\ _Bq_ 

QJ I 

e* 

V 2 

16 q 

/ H 


3B-* 

x> — 2 

3B 2 


— J3 _ 


Ba 

2 

Ba 

3B; 1 


BBu 

2 

= - Br 2 

2 

3e 

2 

Be 

3B 5 > 

H 

= " B '2 2 

BR, 


_ 9 _ 

da 



V 



9B , 

3b, 

- b, b: 2 — - 

Ba 

1 2 3a 

db i 

, 3b 2 


— b . b“ 2 

Be 

1 2 3e 




b 2 2 


3b 2 



Bp 

Ba 



db 2 

Be 


_Bp 

Be 


*-n 0 \p ) *v 0 
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where P n (x) is the Legendre polynomial with argument x of degree n f and P„ (x) is the deriva- 
tive of the Legendre polynomial with respect to the argument. The definition of R is as given 
previously, viz., R n (*) = x n P n ( 1/x). All infinite series are computed by an iterative method, 
with computation of terms ceasing when the absolute value of the ratio of successive terms 
minus unity is less than or equal to some pre-selected tolerance. 


p-‘ H 



- 7] <n “ 2 *(^)” [<1 - r>, ( 7 ) P„_, [<1 - 
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3A, 

3a 




3A, 

Be 


= Aj p 


dp- 


Be 


■ Aj e (1 ~ e 2 ) -1 


+ (1 - e 2 ) 1/2 p- 1 






n = 1 


[(1 - e 2 ) 1/2 ] n _1 



[(1 - e 2 )' 




^2 

Ba 



BAj 

Be 


A . a ™2 

" " 2 


[(1 - e 2 ) l/2 ] 


R„ [(1 - e 2 ) 1/2 ] 


p; [(1 - e 2 ) -1/2 ] 



[(1 -e 2 ) 1/2 J 
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— =-A„ P -‘ -b,p-. |E 


3a 


3a 


*lr- 2 p-’(3b;-b;,U + 2p-. ( 3 b, ^-b/i) 


2 2 


“2 I 


3b, 


-i [3 b, 

Bbi 

\ 1 

3a 

iE -b. 

Bb, 

3a 2 

3a 

■if)} 



2b, 


- A 

Be ” 21 


2 W + e -i 1 3p 


-e (1 - e 2 )" 1 + e" 1 - p 


3e 


P-.fb.p- |f *^-2p- (3b; -b», |f 


"2 2 Be 


t2p ‘‘ (“‘ S 1 - 11 - If * I 1 ’;-- 1 ( lb 7 


Bo Bb. Bb, 

1 - b„ — - - 2 b. — - 


e 3bib 2 p 


2 F Be “ 2 Be 


Be 


.3,3 -3 

+ 2 b 2 P 


(4 + 3eJ) (ir " b2P_1 S) + t b2P ' 2e ( b 2P‘ 1_b i)| 


BA. 


~ = ( i-e 2 )i/2 p- 2 e J_I1 + 2p -i (3b. - b, Z* 

OV, 1 B-n \ 1 3,, 2 2. 


dr >o 


,2\l/2 n -2 J_^l , 


dr >o 


3b. 3b, 


d V 0 2 d V 0 


-|b, P - (. ,lJ) ( b , !5 l . ab , |bf i b » p -» <4 , 3 «■> 


A 22 P* 1 || + J P- 3 e 2 (1 - e 2 ) 1 /2 | - 


= -A .-l^P +1 

3a 


-P - 1 (3b 2 - b 2 ) 
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— p- 2 bj b 2 p- 3 b* (6 + e 2 ) 


Bp 


3b, 


( 3b i -^P-’ b 2 ) 


9 p" 1 b i b 2 + — p- 2 b 3 (6 + e 2 ) - b 2 


Bb, 

3a 


} 
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where 


where 


where 


BA 


22 


-A, 


Be 22 


[ p -‘|£ 4 ,<l 


+ -|-P' 3 e 1 (1 (3bJ -b 3 ) +^P' 2 b, 

-|p-*b< ( 6 *e")j|£, (ib.-fp-bj) ^ 


-]&} 




- 9 p~* bj b 2 + ~p ~ 2 b 2 (6 + e 2 ) 


+~p " 3 e (1 - e 2 ) 1/2 


[ 3 bJ - b§ - 9 p‘» b, b 2 + |p ' 2 b< (3 + e 2 ) J 


BA 
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9B 3 = 3t ?2 

B a B a 
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(1 - V' 2 2 )~ 3/2 V 2 3 + 2 m y m 

m-2 _J 

y m has been given above in the section titled "Mutual Constants. 

CD 
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m-2 
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(1 
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1 ! > 2 2 2n (n ! ) 2 '° 
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= - 3 », p-' || * <1 - .«)"■ p-> ^ K 1 -•’)■'*] -jf 

n-1 

3D n /3a is computed as follows: 
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If n is an even integer, then 
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(- 1 )*-■ (i -m) 





2 m 


1 2 m 





If n is an odd integer, then 


Ba 
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2 i + 1 


Ba 
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^3 
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B e 
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[(l-e 2 )' 1/2 ]-^] (n + 2)[(l- e2 ) 1/2 ] n + 1 

n = 0 


u n P nt2 [d-< 


where D n has been given above in the section titled "’Mutual Constants,’" and where BD n /Be is 
computed as follows: 


If n is an even integer, then 
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If n is an odd integer, then 


^£n 

Be 


BD 


21+1 = - 2 op" 2 Bp 


Be 
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l) 1 -" (i -m) 


m “ 0 


(§T , - , -£) 


2m + 1 
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m = 0 


^m + l) 





2m+ 1 


2m + l 


BA- r — i 

I7 = d- e2 ) 1/2 P' 3 2_, R "^ [d- e2 > 1/2 ] 

° n = l 

where 3 D o /3t? 0 is computed as follows: 

If n is an even integer, then 




3D 3D,. 3 b, 

= 2 p -»_2 V (- 1 ) 

3t?0 37,0 3t) o 




If n is an odd integer, then 
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3D 3D 


Bl 7 0 Bt >o 


(- l) s ~ m (2m + 1) ^ 2<i ' n> 
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>-l.|P + (l-e*)l/2 p-« 
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+ 2p" 2 (4 + 3e 2 ) 
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b 2 P“ 
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~ P 
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+ (1 - e 2 ) 1/2 p 4 e 




33 


2 p~ 2 (4 + 3 e 2 ) ^i-b^+c 2 ^ 


1 ( 

3 + ie 2 \ 
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3 + 4 / 
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B A , 
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BA 
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3 3b i 1 /, 1 2 \ 3b 2 \ 
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BA 
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Be 
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'd7] c 
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BA 


Ba 


33 = -3A„ p" 1 ^ + (1 - e2 ) 1/2 P" 4 e3 
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1 Bb l 1 , Bb 2 

12 Ba ‘ 3 P 2 B a 


Bp 

a a 


[§ p-2 (j b2 +c2 ) ~T2 p ‘ lh ] 


BA 
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Be 


= - 3A 33 P' 1 ^7 + p " 

4 e 2 

[3 (1 

+ (1 - e 2 ) 1/2 p~ 4 e 3 < 

f 3 p 

[’p- 

l de 

L 3 
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BA 
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BA 
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Be 
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bb : 1 . 
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~J7 = 2^ V '° l (a +b » +A 1 +c2 ^0 A 2 B 1 B 2 1)_1 (( a 2 - a 3) 1/2 

I B a 

/ Ba Ba,\ 

+ ( a 2 ~ a j) - 1/2 A 2 Bj 1 ( a 2T f- a 3 
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17 Be ^o® 1 ® 2 Be 
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The following time- independent partial derivatives are computed only if the differential 
correction is carried through terms of second order: 


= - 3 A n p' 1 ~|(1 - e 2 ) l/2 p' 3 e b 
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. _ Bb, Bb 
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BA 


23 


3a - a =3 p_, lf b = 
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Bb 
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BA *» - A 
3a 24 
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TT = *» ( 4 b i‘^T" 5 If) * 5§6 b ' s b ; ** [ 4 - <* - eI >' ,/t e ‘] 


BA^ 

3t )o 


* ^24 ^*2 
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H 


DIFFERENTIAL CORRECTION: TIME-VARYING PARTIAL DERIVATIVES 
WITH RESPECT TO ENERGY-MOMENTA VARIABLES 

We now compute the following partial derivatives of time- dependent parameters with re- 
spect to the orbital elements a,e > and v 0 . We shall later compute the partial derivatives of 
these same parameters with respect to the remaining orbital elements P l , P 2 , and P 3 . 
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(1 - e 2 ) sin E (sin v )* 1 (1 - ecosE )" 2 -EE 

Ba 


3 M S Bv 0 By, 
Ba + Ba + Ba 


f I /BM Bv. Bv,\ 

[(1 - e 2 ) sin E £ + sin 2 E J (sin v )' 2 (1 - e cos E )' 2 - 


= (1 - e 2 ) sin E (sin v )~ 2 (1 - e cos E ) -2 


BE 


BM s Bv 0 Bv, 

.H + H H 


BM s 

B v ° 

+ !Ei 


Ba 

Ba 

Ba 

Ba 

BM s 

3v o 

Bv, 

Bv, 

Be 

+ Be 

+ Be 

+ Be 

BM s 



+ !El 

? v 0 

d % 

+ 

3t ) 0 


' = j^2 ‘ ^ ^ B 2 1 [^1 C ° S < 2 ^s + + |" ^ Sin ( 2 ^s + 2 ^o) 

’ 0) ]} “ \ S 1 +(a2 - a 3 )_1 


— q sin ( 4 ^ s + 40 o ) 


Ba Ba 

+ (a* - a*)* 1 a 


II _ b : 1 ^ 


(- 2 a ,)' 2/2 ( a 2 - a 2 ) 1 2 t ;' 2 B * 2 


Bv, 


BA. Bv 


A2 TT + " 2 TT + bT (A22 cos v ' + 2 A22 cos 2v<) 


- (A 2 i v j sin v' + 4A 22 v 1 sin 2v' 


3 A 23 cos 3 v* 


4 A 24 cos 4 v' ) 



4 \\ COS V # 



Ba 


1 A. . B A 2 -j , ' ^24 

4 2v. cos 2v' + sin 3v # — 4 - sin 4v — — 

da Ba Ba 
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7 q B >' (if - 7 qB '.‘ 5)] [*■ ”* < 1 2 *. * 2 *.> * f 


q 2 sin (2^ s - 2y 0 ) 


— q 2 s in(4i/> s + 40„) 


if r 

- q 2 BJ 1 <cos (2</> s + 2i/> 0 ) -i- - |2^ sin (2^ s + 2^) 


if •'* cos (2^ + 2</' 0 ) + y| ‘l 2 cos ( 4 ^ + 4^ 0 )j ^ ^ 


3 

+ _q 


1 


sin (2t/> s + 2i/> 0 ) - - sin (4</* s + 4</< 0 ) 


•>] £} 


B02 

Be“ 


B 2* 


cos (2'/' s + 2^ 0 ) + i q 2 sin (2<p s + 2\p 0 ) 


3 

— q 2 sin (40, + 40 o ) 


>]} 7 * <“• - (“* S - “» §) - B V 5] 


(- 2aj)- 1/2 (a 2 - a 2 ) 1/2 T^Bj 1 


Bv, 


3A 2 BV 1 


*2 sT + V 2 17 + 17 (Aj, COS v' + 2 Aj 2 cos 2v*) 


-(A 21 v, sin v' + 

^ ^22 Vj sin 2v' - 3^ cos 3 v' - 4 A 24 cos 4 v') 


/BM, 3vA 
\Be Be/ 


, aA 2! 


BA, 


BA, 


+ V cos v' — — + 2v cos 2v' 11 + sin 3v' — H + sin 4v # — — 


. 3A 24 


Be 


Be 


Be 


Be 


1 q B ' 1 [fl - I q B * 1 — \ 

2 2 \3e 2 2 Be / 


0 X cos (20, + 2i/> 0 ) + .1 q 2 sin (20, + 20 o ) 


2 

q 2 sin (40, + 4^ 0 ) 


+ I q 2 Bj 1 <! cos (20, + 20 o ) 


B^ 

Be 


-- 20 x 


sin (20, + 2 0 O ) 


•3 3 / B 0 B 0 A 

| q 2 cos (2 0, + 20„) + — q 2 cos (40, + 40 O )J + — J 


3 

+ 4 q 


sin (20, + 20 o ) - — sin (40, + 40 o ) 


4 >] &} 
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B0 



3^ 

CM I 

| 

- 

Ba 

Ba 

Ba 

3a 

3a 


_ 


3^ 

+ rr — 

Be 

Be 

Be 

3e 

Be 

Bi p 

30 s 

#0 

+ 

3^ 

3^2 

~o 



+ 3% 

3r, 0 


This completes the computation of the partial derivatives of the uniformizing variables E,v , 
and ^ with respect to the orbital elements a, e, and % when the calculation is followed through 
terms of the second order. If, however, second-order precision is not necessary, we can 
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eliminate the terms with the subscript ”2" (thus omitting all partial derivatives of M 2 , e 2 , v 2 , 
and ' 4 > 2 ) , and the above partial derivatives of the uniformizing variables reduce to the following: 

BE = + 

3a 3a 3a 


3E _ BE 

Be Be Be 


BE _ B8 BEj 

^ o H ^ o 


3v 

Ba 


BM s ^ Bv 0 ^ Bv t 
3a 3a 3a 


3v = dv o Bv i 

Be Be Be + Be 

3v _ 

*V 0 Bt, 0 

Bi p _ ^ 

3a 3a 3a 3a 

3 ± _ 

Be Be ^ Be * Be 

= ifi 3 ^> B ^i 

Bn 0 Bt5 0 B7, o 

We now continue with the necessary equations preparatory to the partial derivatives of the 
orthogonal co-ordinates X , Y , and z. 

^ = (1 - ^) sin i// (sin y) _1 (1 - 7)2 sin 2 v //) -372 


4^- = (1 - rfi) sin i/; (sin y) _1 (1 - t) 2 sin 2 i/>) -3 2 

-?L= (1 - ’ij) sin Wsin y)* 1 (1 - ^sin 2 i|i)' 3/2 3^- 

ot? 0 cn 0 


- t ) 0 cosi/i sin 2 >/■ (sin y )' 1 (1 - t ; 2 sin 2 i/') 3/2 
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30 

3a 


a - ngr 1 " a * + b,,m 3 .2,-4 


^[(“a " a i 


; 3 V + 32 T5 ’S' s in 20 ||(a| -a 2 )*" 2 „ 

3a 


”• (“■ 17 - “» ir)] * “»<“5 - “D- /1 % [<• - *>" 


(1 is. 

3a 


- V a - % 2 r 1/2 (i - ^ 2 r 3/2 v ^ + b 3 |i + * 

3a 3 3a 3a 


3 O c *J 

^2 Sln 2^ ^ + ^ 


r B0 1 

3a J 


- c 2 


bT + Te 1,0 7,-2 cos 

l < - 2 *> ) " / ’ IT * “> <- 2 “. ) ' S/2 5]<V * *1" v , 

+ Ajj sin 3v + A 34 s i n 4 v) - c 2 a 3 ( - 2a t )' 1/2 £ 

+ 3 A 33 cos 3 v + 4A 34 cos 4v) |l- + sin v + sin 2y ^32 

3a ^ 


A 32 sin 2 v 


( a 3 + A 31 cos v + 2A 32 cos 2v 


3a 


3a 


+ sin 3 v 


3A 33 BA .4 3A,1 

+ sin 4 v —11 + v -1 
3a 3a J 


Be (1 " V o yU2 (1 " ’’j 2 )' 172 V + B 3 0 vl V~ 2 4 sin 2 <P 


(a 2 - a 2 )' 1/2 tj 0 ^L 


a 3 ( a 2 - a2 )‘ 3/2 


/ Ba Ba \ 

M* 2 ^ “» J?) 


+ a 3 (a 2 - a 2 )" 72 


(1 - T ) 2 )-!/ 2 (t _ 7 ^“ 2 )~ 1/2 i>C 

3e 


- ^)' 1/2 (1 - ^ 2 )- 3/2 y ^ + B^ + 4 , 

3e 3 3e 3e 


3B, 


•3 2 — cr 

Tjf 77 5 




? o 7 7 2 5 sin + ^ ^0 77 2 4 cos 20 


3e 16 'o ; 2 


30 

3e 




<-2 V — ^ (A; 


> 3 v + A 31 sinv + A 32 sin 2 v 


+ A 33 sin 3v + A 34 sin 4v) - c 2 a 3 (-2^)- 1/2 J^(A 3 + A 31 


cos v + 2A 32 cos 2v 
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B<£ 


+ 3A 33 cos 3v + 4A 34 cos 4v) + sinv + sin 2v^^- + sin3v ^-+sin4v 

° e 3 e Be Be Be Be . 

= [(l-^r 1/2 (l-^V 1/2 y + B^ + ^^^sin2^(a*-a2)- 1/2 ^ 0 
- «*,<«* - -jr" 2 % (a 2 - a 3 ^ + «,(«» - a’)-"’] 

!>- 1/2 ^[(i -ir ,/2 <i - ^ 2 r 1/2 |3L--- 3 a -^>-» /2 (i _ 7?2 - 2 )-3/2 y i^ + B +v ,25l 

L d7, o 3t ? 0 H 37 7 0 


- (r.2 _ „2>* 

~3 V ~2 ~3' 


Bp. 


2 3 


■f ^ V~ 2 5 sin 20 +Jj-vl V' 2 4 cos 20 + 7*> (1 - ^)" 3/2 (1 - V 2 2 )' l/2 y 


+ lT r, o 77 2 4 sin 20 


r Ba, Ba, 1 

C 2 j^(- 2 a,)* — i + a 3 (-2V- J/1 ^j(V^. sinv 

+ A 3 2 s in 2 V + A 33 sin 3v + sin 4v) - c 2 a 3 ( - 2a t ) 1/2 ^(A 3 + A^ cos v + 2 A 32 cos 2v 


Bv BA—.* BA^i% BA* ^ 4 

+ 3 A,- cos 3 v + 4 A~ 4 cos 4v) — — + sin v — — + sin 2 v — — + s i n 3v — — ~ + s in 4v - + v 

'*** OT) ^ /4*n H-ti 7^-n 


3t *> 


3 ' % 




£■ x [< 


B 0 l_ Y 30 
Ba J Ba 


(p 2 + c 2 )~ 1 p|aesinE + 1 - e cos - (1 - p 2 s i n 2 i/>)~ 1 p 2 sin 0cos ^ 

~ X J^(p 2 + c 2 )" 1 pa|e s in E - cos E^ - (1 - p 2 sin 2 p )~ x p^ sin ^ cos p - 

-7 4 


B<£ 

17 


(p 2 + c 2 )" 1 p a e sin E 


BE 


- (1 - Vg sin 2 0T 1 T) 0 sin 0 ^tj 0 cos 0 + sin 0 


Bl// 




B<p 


BY 

Ba 


Y j(P 2 + c 2 )~ x p ^a e sin E + 1 - e cos 

- (1 - r, 2 sin 2 0)- 1 tj 2 sin 0cos 0 |^j + 


X 


30 

Ba 


— = Y O 2 + c 2 )- 1 
Be L 


a e sin E 


BE 

Be 


cos 


E ) 


- (1 - ^ sin 2 0)* 1 


V 2 sin 0 cos 0 


B0 

17 




X 


B0 

Be 
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= Y (p 2 + c 2 )* 1 pae sin E 

I Bt »o 


- (1 - ^ sin 2 i/-)- 1 7) 0 sin \P (ri 0 cos 4> + sin A 1 + X 


^- = p 0 (1 - e cos E) f sin >/< + a cos ip ^ j + a e -q 0 sin 4> sin E 


• = ri Q (1 - e cos E) a cos \p + a tj 0 sin *p ( e sin E ^ — cos E 


— = a(l - e cos E) (sin ^ + V 0 cos r^~ j + a e p 0 sin <// sin E — 


DIFFERENTIAL CORRECTION: TIME-VARYING PARTIAL DERIVATIVES 
WITH RESPECT TO ANGLE -EPOCH VARIABLES 

We now compute partial derivatives of the time-dependent parameters from the orbit 
generator with respect to the orbital elements / 3 ,, A,, and A 3 . This procedure is analogous to 
the one followed in the preceding section. Whenever a partial derivative with respect to P 3 is 
not given, it is assumed to be zero. 


= ‘ 27Tl 'l c2? ?0 “j' B 1 B 2 l 


3/3 t 2vv i 


30 

— - = 2 77 V a’ 1 A" 1 (a + b + A x ) 

3 / 3 , 


— - = (1 - e' cos g) 1 

3 / 5 , Wx 


=(!-•' cose ) _1 


/ "iO 

— = (1 - e 2 ) sin g (sin v')* 1 (1 - e cos g) -2 — - 
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38 3M S 


(1 - e 2 ) sin 8 (sin v')' 1 (1 - e cos g) -2 — — - - 


3/3, 3/3, 




^-=(-2a irl / 2 («g -«*)!/» ^ 


oM, Bv n 

— =-( 8 +^)-* jCAj + c 2 t£ AjBjBJ 1 ) — 


1 ( 0X P 

•— c 2 t^ (- 2a J ) 1/ ' 2 (a 2 - aj)- 1/2 cos(2^ s + 2^ 0 ) ( — ^ + 


3^ s B0 O ’ 


3/3, 3/Sj, 


3M, 

-_=-(a+b 1 )- 1 (A 1 + c 2 t^A 2 


B i B , ) 


, /B0 3i// 

— 2 c2 ^ (- 2a i) 1/2 ( a 2 - a 3>' 1/2 cos < 2 ^ + ^o) \^ + 


BE, 3M, 

— — = (1 - e' cosS)' 1 — — [l - M, e ' ( 1 - e' cos g) 2 sin8] 


M. e' (1 - e' cos g)~ 2 — - [sin g 


+ — Mj (1 - e' cos 8)' 1 cos 8 - — M,e' (1 - e' cos 8)“ 2 sin 2 8] 


3E, , 3M 

= (1 — e f cos g) 1 — — 

3/S 2 1 ; 3/3, 


- = (1 - e* cos g)- 1 — - [1 - M, e'(l - e' cos 8) 2 sin 8] 

< v ' ~A R 


-M,e'(l-e' cosg)' 2 [s in 8 + M,(l - e' cos 8) 1 cos 8 

°A> L 2 


Mj e' (1 - e' cos 8)* 2 sin 2 8^ 
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3 v 

(1 -e 2 ) sin (g+E,)(sin v*)" 


1 


[1 - e cos (g + E ,)] -2 


' 3 g BE/ 

M + 


' 3M s 3 v o\ 

+ *A/ 


— L- (1 - e 2 ) sin(g + E.) (sin v ")" 1 [1 - e cos(g tEj )]- 2 

B/ 3 2 


Bg_ 

B/3 2 + B/3 2 


'BM, BvA 

^ + B/y 


‘dxp 

- 11 = (- 2 a ,)-" 2 
1 


(a 2 _ a 2 )l /2 T)" 1 B" 1 



+ (A 21 cos v # + 2 A 22 cos 2 v') 


BM S BvA" 

M + 


+ 1 q 2 B ' 1 cos (2 *p s + 2 i/> 0 ) 


M + B/ 3 j 


B'/'j 


1/2 


( a i 


a2) 1/2 77" 1 B“ 



+ (A 21 cos v' 


+ 2 A 22 cos 2v # ) 


BM S + Bv q \ ~ 

+ B/ 3 s A 


+ j q 2 B * 1 cos ( 20 s + 20 o ) 


B'/'j B 0 o \ 
B/3 2 + B/3 2 J 


The following time- dependent partial derivatives with respect to the orbital elements p l , 
p , and p are computed only if the differential correction is carried through terms of second 
order: 


bm 2 

3 A 


- (a + bj )* 1 



+ (A u cos v' + 2 A 12 cos 2v') 


(l ^ jM 
Wl + B^J 


+ c 2 rj* (_2aj)l/2 (a 2 _ a 2)-l/2 


B’/'j 


Bi Wi " i cos {2Kps + 2 ^ o) ^ 


B’/'j 


+ </>! sin (20 s + 20 o ) 


'B’/'j + Bi// 0 \ 

B/ 3 , + B^/ 


-|q 2 cos (2 0 s + 2^ 0 ) 


B’/'j 30 o \ 

M + 30 J 


+ 


16 


q 2 cos (4^ s + 4 </> 0 ) 
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3M 2 

^ f =-(« + b , r 1 

'"'^'2 



3V 1 

TT- + (A,, cos v' + 2Aj 2 cos 2v') 



+ c 2 ^ (- 2a 1 ) 1/2 (a 2 - a 2 )* 1/2 js , ^ - y cos (2 0 S 4 20„) 

+ ^ sin (20 g + 20 o ) ^—4 ^g-j - -q ! cos (20, 4 20 o ) ^ 4 


**o\ 


+ 


16 


q 2 cos (4 0 s 4 40 o ) 



3E, 


[1 


e* cos (8 4 Ej)] -1 



- Mj [1 - e' cos (8 4 E,)]* 2 


e' sin (g 4 Ej) 



I! 2 

^2 


[1 - e' cos (8 4 Ej)]* 1 


w 2 


-M 2 [1- e'cos(8 4 Ej ) ] * 2 e' 


_3E_ _ 38 

3p, _ 3/3j + 3/5, + 3/5, 
3E _ 3c ^®2 

~ 2 ~ ~ 2 + W 2 


sin (8 4 Ej) 


'38 

W 2 


3EA 

3 / S 2 / 


— — = (1 - e 2 ) sin E(sin v)* 1 (1 - e cos E)* 2 

B6, 


— - = (1 - e 2 ) sin E(sin v)" 1 (1 - e cos E) 2 
^2 


3E 

/3M S 

3v o 

+ ZM 

3/3j _ 

V01 

f 3P, 

BPjJ 

3E 

/3M S 

3v o 

3vj' 

3/3 2 “ 

W 2 

+ ^2 

+ ^ 



3v 

BM s Bv Q 3 v 2 


3^j 

3Pj ' 3/5, ' 3/5, 3.5, 

- 

OJ 

< 

3M S 3v q 3vj 3v 2 


1 

1 M 

3/S 2 3p 2 + 3,6 2 + 3/5, 
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Ip (-2 - X )- 1/J <«» -«>)»/* 


Bv„ 




3v 

+ -^-r- (A 21 cos v* + 2A 22 cos 2 v') 


- (A 21 v, sin v' + 4 A 22 Vj sin 2v' 


3 Aj 3 cos 3v* - 4 Aj 4 cos 4v') 


3M s 9v oV 

M 3/3,/ 



,cos ( 20 s + 2«/> 0 ) 


3^ 


2'p 1 sin (2> p s + 2i/< 0 ) 


-■j q 2 cos (2^ + 2^) + — q 2 cos (4 <//, + 4</< 0 ) 

4 16 



Bvi/ 


( a 


a 2 )l /2 




I 1 A, 


Bv„ 


Bv, 


3/8, 3/8, 


(Ajj cos v' + 2A 22 cos 2v') 


- (A 21 Vj sin v' +4A 22 Vj sin 2v' _ 3 A J3 cos 3v' 


4 A 24 cos 4 v') 




cos (2</< t + 2 </< 0 > 


3^1 

3/3 2 


2^ sin (2 v// s + 2^ 0 ) 


3^ 

4 


cos (2V^ S + 2 0 O ) + 


_3_ 

16 


q 2 cos (4^ s + 4^ 0 ) 



B/\ + B/3 X + B/3 X + B^ x 


Bi/i _ B^ 0 B^ B^2 

B^~1^ + B^ + 3/3 2 + B^ 2 

This completes the computation of the partial derivatives of the uniformizing variables 
E , v , and \p with respect to the orbital elements j3 i , fi 2 , and /3 3 (those with respect to p z are all 
zero) when the calculation is followed through terms of second order. If, however, second- 
order precision is not necessary, we can eliminate the terms with the subscript "2” (thus 
omitting all partial derivatives of m 2 , E 2 , v 2 , and ), and the above partial derivatives of the 
uniformizing variables reduce to the following: 

BE Bg BE i 
3/3,' 3/3, + 3/3, 
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BE _Bg - E i 

a .£ 2 a /? 2 a /? 2 

_?V__ 3v Q 

a/?, = ^ + a/3j + ^7 

_Bv = ^_ ^t_ 

B/3 2 ' B/3 2 + B/3 2 + 3/3 2 


30 ___^ ^ 2 _ 

a/?j " a/s, + a/3j + b^ 


By _ ^b_ _^L 

W 2 ~ W 2 + a ^7 + a^ 

We now continue with the necessary equations preparatory to the partial derivatives of 
the orthogonal co-ordinates x, Y, and z. 


= (1 “ ^o) sin ^ ( sin x)' 1 (1 - ^ sin 2 0)' 3/2 

d Pi 3 / 3 j 

T^= (1 - sin 0 ( si n v)" 1 (1 - V* sin 2 0)“ 3/2 

'2 d /°2 


By 

u Pj 


; ( a: 


a ?)~ 1 


(1 - 


rtr 1 ' 2 


a 


- 2 \ - 1 /2 


ay 

^1 


+ B, 


^-+2- ri 
B/3j 16 0 


cos 20 



- c 2 a 3 (« 2 a J )" 1/2 (A 3 + A 31 cos v + 2 A J2 cos 2 v 


+ 3A 33 cos 3v + 4A 34 cos 4v) 


Bv 



,( tt 2 


* 3 >‘ 




(1 


t>* 2 \— 1 /2 


ay 

a/s 2 


d 0 

a^ 



^ r ,* 4 cos 2 ^ 
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- c 2 a 3 (-2a 1 )” 1/2 (A 3 + A 31 cos v + 2 A^ cos 2v 


+ 3A 33 cos 3v + 4 A^ cos 4v) 


Bv 

w 2 


_ . 

- (i - ^ sin2 'Z')' 1 ^ sin 'P cos 'Z’ ||^j - Y ||- 

^X^cVpaesinE || 

- (1 - sin 2 i/-)" 1 7) 2 sin \p cos «/> BBj - Y 


3X _ Y 

3 A'“ 



• „ be 

p a e sin E — - 

B/^j 


- (1 - V„ sin 2 i//)' 1 7] 2 sin 0 cos <// + X -|^- 


- =Y [ ( 
^2 L 


Be 

(p 2 + c 2 ) - * 1 p a e sin E — - 
^2 


- (1 - tj 2 sin 2 0)-‘ 7) 2 sin 4> cos 0 ||^j + X 


B0 I l Y B0 


BY 

B& 


= X 


bz r 

*r , n < 


(1 - e cos E) cos 0 + e sin 0 sin E 

B^ o/3, 


f] 
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3Z 

^2 


a [j 1 


/ B0 

e cos E) cos 0 

BA 


+ e sin f sin E 



BZ 

^3 


0 


THE EQUATIONS OF CONDITION 

Now that we have found 3x/Bq , 3Y/dq. , and 3z/3q. for the Izsak orbital elements q. (i = 1, 
2, . . . , 6), we can complete the differential' correction process by deter mining the equations of 
condition. First we expand and substitute into the matrix relation given in the section titled 
"Analytical Procedure of Differential Correction". The matrix relation, when expanded explicitly, 
yields the following eighteen equations: 


dX u 


— — - cos 
da 


0 

“x 


BX 

da 


+ sin 0 x 


BY 

Ba 


Ba 


sin 0 sin di 


BX 

■=r + COS V sin 

Ba * 


Ba 


+ cos 6L 


Ba 


— - sin 0 x cos d D - -COS 0 x cos 6 , D — + sin 6> D _ 


Ba 


Ba 


Ba 


Be 


, BX . 

cos 0 — - + sin w 
Be x 


BY 

Be 


3Y m ... . . BX BY 

- — - - sin r - sin : + cos y sin _ - — 

d e * u d e * u d e 


BZ 


ni , ^ BX , BY . BZ 

■ : sin w cos t - cos 0 cos — 4 sin 

Be D Be * x D Be D Be 


,X M , 3X dY 

cos vs. — — + sin y; 


37 / n 






■ j * /-> BX ; - O’ I ^ 

“ sin 0 x sin d D — + cos 0 x sin $ D — ~ + cos d D — 


BY 


BZ 


Bn 


3t)„ 


37 J n 


A = sin cos * D |L - cos ^ cos * D |L + sin |L 


i BX BY 

_ = cos _ + sin — 
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3Y m 

- sin 0 s i n 0. 

3/3. 1 


3X 


+ cos 4> x sin # D 


3Y a 

w t + cos D 


3Z 

3/3. 


BZ m 

— = sin </> x cos 


3X 

3/3. 


- COS l/> x COS 


by 

3/3, 


+ sin (9 d 


BZ 

3/3, 


M , ■ i 

= cos w + sin \p 

3/3 2 x 3/3 2 


3Y 

3/3 2 


J^u 

3/3 2 


= - sin \p x sin 6 d 


BX 


+ cos if x S in # D 


BY 

3/3 2 


+ cos 0 D 


BZ 

3/3 2 


BZ m 

^7 


sin 0 x cos 


BX 

3/3 2 


— COS 0 cos 

r x D 


BY 

^2 


+ sin 0 D 


BZ 


= COS 0 

3/3 3 


3X 

3/3 3 


+ sin i// x 


BY 

^3 


3Y., 


^3 


s i n 0 x s i n B d 


BX 

^3 


+ cos sin 0 B 


BY 

^3 


M • i n UA , O I 

— — = sin ii/ v cos tL — — - cos (p cos 6L 

3^3 x D 3 fi 3 x D 3/3j 

The last two equations have only two terms on the right-hand side because of the fact that 
3Z/3 /5 3 = 0. We can now write out explicitly the twelve coefficients to be inserted into the equa- 
tions of condition: 


+ h. *0* + .!b 

3a ” dXy 3a 3Y M 3a ; 3 Z Mi a a 


3L. 

3e = 3X„ 3e + 3Y„ 3e + 3Z M 3e 

34 _ 31^ 3X„ + 3^ 3Y„ t l3L c '3 Z m 

3% = aX M ^ r 'o + 3Y m B77 o + 3Z m|H 


3L 


3 


3L c 

3X m + 

3L c 

3V m 

3Lc 3Z m 

?X m 

+ 

3V m 

3P, ’ 

' 3Z m 3,3, 


3Lc 3L e 3 X m | 3L C 3Y M | 3L e 3Z M 

3/3 2 = 3X„ 3/3 2 + 3Y m 3p 2 + 3Z„ 3,6 2 
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« 


3L C 

3L c 


BL c 

3Y„ 

BL C 

BZ„ 

^3 

9X„ 

^3 

+ 9Y m 

9^3 

+ BZy 

9*3 

9M C 

3M C 

bx m 

BM 

. £ 


3M_ 

C 

9Z* 

B a 


Ba 

9Y m 

Ba 

+ bz; 

Ba 

BM C 

9M C 

9X m 

3M_ 

C 

3Y « 

BM C 

BZ„ 


oM c oM c 3 X m 3M c BY„ 3M c BZ„ 
~ 0 ' BX^Bt^ + + 9t7 


BM C BM c BX m BM c 3Y m 3M c BZ m 


B/3, ' 

' 9X m 

B/3, 

+ 9Y M 

9/3» 

+ 9Z„ 

B/3, 

3M c 

bm c 

bx m 

bm c 

9 y h 

BM C 

BZ || 

B/3 2 
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Finally, the two equations of condition corresponding to each observation are given ex- 
plicitly by the following: 


AM = M 0 -M c 


BL BL. 

C . C 

A e 

9L c 

T — — 3 -j. — 

c* a c e 

+ 9^o 

dM 0 3M r 

1 + 1 

A e 

BM C 

3 a 3e 


+ 9 Tn 


BL. BL BL 

At *0 + — A A + T/T A ^ 2 + ~ A/3 3 


B/3, 

3M 


BM C 

9/3, 


9/S, 


BM c 

3/8, 


FITTING BY METHOD OF LEAST SQUARES 

We have accumulated a set of 2n linear simultaneous equations in six "unknowns," as 
follows: 


(A L) ; = (L 0 ). - (L c ), ~- 


(AM), = (M 0 ). - (M c )j = 




y 


A< *i 


i = 1 , 2 , ' • • ,n 


where (j = 1, 2, . . . , 6) = a, e , 77 0 ,/3, ,p 2 , p 3 . We regard the Aq ; as "unknowns,” and the 
number n of observations in the set is fixed in advance (see above under the section titled, 
"The Standard Deviation of Fit"). The above equations, written in matrix form, become: 
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[A qj ] - [AL.] 
[Aq } ] = [AM.] 


where the matrices of partial derivatives have n rows and six columns, the matrices of un- 
knowns have six rows and one column, and the matrices of observational residuals have n 
rows and one column. Recall that, in the general case, the observed direction cosines (L 0 ) l 
and <$A 0 ) i have associated weighting factors (w L ) i and (w M ). respectively (see above under the 
section titled, "Computation of Direction Cosines"). 


For purposes of this section, it is unnecessary to distinguish between direction cosines 
L andM or between weighting factors w L and w M . Further, it is not significant, for the present 
purpose, that the constant coefficients in the linear simultaneous equations have the form of 
partial derivatives. In order to simplify the notation in what follows, we combine the two 
matrix equations, each coefficient matrix having n rows, into a single matrix equation where 
the coefficient matrix has m = 2n rows. Then the matrix of constant terms (i.e., observational 
residuals) also has m rows. We rewrite the above two equations in the simple general form: 


AX = B 


where A = [a..] has m rows and six columns and represents the coefficient matrix of partial 
derivatives, x = [x.] has six rows and one column and represents the matrix of unknowns, and 
B = [b.] has m rows and one column and represents the matrix of observational residuals. 

The number m of equations we obtain by expanding the matrix relation is generally much 
greater than the number (six) of unknowns, and since the observations contain inherent random 
and possibly systematic errors, no exact solution of the simultaneous set exists. According to 
the principle of least squares, the values of the unknowns x . which are preferred are those 
which cause the sum of the squares of the residuals after the fit to be a minimum. The so- 
called "residuals after the fit" are calculated by substituting the approximate solution for the 
xj in the matrix x, and subtracting the matrix AX from B. When the equations of condition 
have different weights, the least-squares solution is that which minimizes the sum of the 
weighted squares of the residuals after the fit, where each square is multiplied by its corre- 
sponding weight. 

The least-squares criterion is satisfied by reducing the m equations of condition to six 
equations known as normal equations. This procedure is performed as follows, in which we 
adopt the usual notation for matrix elements: the first subscript denoting the row number and 
the second subscript the column number. The first normal equation is obtained by multiplying 
the first conditional equation by w x a n , the second by w 2 a 21> the third by w 3 a 31 , etc. and sum- 
ming the resulting m equations. The second normal equation is obtained by multiplying the 
first conditional equation by w x a 12 , the second by w 2 a 22 , the third by w 3 a 32j etc. and summing 
the resulting m equations. If we repeat this process six times, we obtain the six normal equa- 
tions. It is seen that this process is equivalent to pre- multiplying the matrix equation AX = B 
by the weighted transpose of the matrix A, where the rows of the transpose are multiplied by 
the corresponding weighting factors. The set of normal equations can be represented by the 
new matrix relation CX = D , where 


a ki a kj w k (*’ j = i ' 2 ' >6) 


m 
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and 


d i ° ^ ^ b k W k (‘ = 1 - 2 -' • •' 6 > 
k = 1 

Of course, if the weighting factors are not present, i.e., w k = 1 for all k, the elements c. 
are precisely those of A T A and the elements d. are precisely those of a t B. Here the super- J 
script "T" indicates the transposed matrix. 

We now have a system of six equations in six unknowns, since C is a square (and sym- 
metric) matrix. In order to solve this system, we use a method known variously as the Gaus- 
sian elimination method or the method of pivotal condensation. This has the effect of reducing 
the square matrix to an upper triangular matrix (i.e., all elements below the principal diagonal 
are zero) which represents the same solution for the x . . 

To begin this process, we choose the element of the first column of matrix C greatest in 
absolute value, say c kl . We then divide all the elements of the k th row (the "pivotal row") by 
the so-called dominant element (or "pivot"), c kl . This done, we exchange the corresponding 
elements of the pivotal row with those of the first row. The leading element c n of the matrix 
is now unity. We now replace all the elements in each row beginning with the second row by 
the following procedure: multiply all elements in the first (pivotal) row by the element in the 
first column of each row successively and subtract this product from the corresponding ele- 
ment of the successive rows. Mathematically, this is indicated by: 

c ij = c ij - c ii c ij (i = 2 ’3,- * ; j = 1,2,* • *,6) 

Since c u = 1, it is obvious from this equation that c u = 0 for all i = 2, 3, . . . , 6. That is, all 
elements in the first column except for the first (diagonal) element are replaced by zeros. 
Essentially, we have added suitable multiples of the pivotal row to all the other rows so that in 
each resulting row the element in the first column vanishes. 

Consider the matrix with five rows and five columns obtained by deleting the first (pivotal) 
row and the first column. Now select as a new pivotal element the largest element in absolute 
value in the new first column of the five-by-five matrix, and repeat the entire process with re- 
spect to the square matrix of order five. 

Continuing in this manner, we have finally a single non-zero (diagonal) element in the last 
row. The procedure is completed by dividing this final row by the diagonal element. The result 
is an upper triangular matrix with ones along the principal diagonal. Note that all operations 
described above to be performed on the original square matrix C are elementary row operations 
(i.e., an operation belonging to one of the three following types: the interchange of any two rows; 
the multiplication of a row by any non-zero constant; the addition of any multiple of one row to 
any other row). Thus, the triangularization process does not change the solution to the simul- 
taneous set of linear equations as long as the operations performed on matrix C are performed 
in an analogous manner on the elements of the column matrix D. This can most readily be done 
by augmenting the six-by-six matrix C by a seventh column composed of the elements of D. In 
practice, the six-by-six matrix C is further augmented by a six-by-six identity matrix placed 
in columns eight through thirteen. The purpose of this is to determine ultimately the inverse of 
the coefficient matrix C, from which we may easily find the standard errors of the least-squares 
solution for the Aq . Note that the various columns of C _1 can be found in succession by solving 
the matrix equationCX = l i for the column matrix X , where I i represents the i th column (1 - 1, 
2, . . . , 6) of the identity matrix of order six. We can thus view the six columns of the identity 
matrix placed in the augmented six -by -thirteen matrix as constant right-hand side column 
matrices replacing B in successive least-squares solutions of the matrix equation. These suc- 
cessive solutions are determined simultaneously in the Gaussian elimination method simply by 
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forming the augmented matrix E and performing the elementary row operations on all thirteen 
columns. The augmented matrix E appears as follows, after the normal equations are deter- 
mined, but before the elementary row operations are begun: 


u 

C 12 

C 13 

C 14 

C 1S 

C 16 

d i 

1 

0 

0 

0 

0 

0 

21 

C 22 

C 23 

C 24 

C 25 

C 26 

d 2 

0 

1 

0 

0 

0 

0 

31 

C 32 

C 33 

C 34 

C 35 

C 36 

d 3 

0 

0 

1 

0 

0 

0 

41 

C 4 2 

C 43 

C 44 

C 45 

C 46 

d 4 

0 

0 

0 

1 

0 

0 

51 

C 52 

C S3 

C 54 

C 55 

C 56 

d 5 

0 

0 

0 

0 

1 

0 

61 

C 62 

C 63 

C 64 

C 65 

C 66 

d 6 

0 

0 

0 

0 

0 

1 


After the triangularization process, the augmented matrix E is transformed to a matrix 
(call it f ) of the following form: 


1 

f 12 

f 13 

f ,4 

f ,5 

f ,6 

g n 

g 12 

g l 3 

g 14 

g l 5 

g 16 

g l 7 

0 

1 

f 23 

f 24 

f 25 

f 26 

g 21 

g 22 

g 23 

g 24 

g 25 

g 26 

g 2 7 

0 

0 

1 

f 34 

f 35 

f 36 

g 31 

g 32 

g 33 

g 34 

g 35 

g 36 

g 3 7 

0 

0 

0 

1 

f 4S 

f 46 

g 41 

g 42 

g 43 

g 44 

g 45 

g 4 6 

g 4 7 

0 

0 

0 

0 

1 

^56 

g 5 1 

g 52 

g 53 

g 54 

g 5 5 

g 56 

g 5 7 

0 

0 

0 

0 

0 

1 

g 6 1 

g 62 

g 63 

g 64 

g 65 

g 66 

g 6 7 


The first six columns of F represent the triangularized coefficient matrix, and the re- 
maining seven columns represent successive constant right-hand side matrices, each of which 
is associated with a particular column solution matrix. At this point, it is only natural that we 
augment the column solution matrix x (corresponding to the seventh column of F only) to a 
six-by-seven solution matrix Y, which contains X as its first column. The remaining columns 
of Y will contain the inverse of the coefficient matrix C of the normal equations. 

We can now write explicitly the set of linear simultaneous equations in the triangularized 

form. 


y ii + ^12 y 2i + ^13 y 3i + ^14 y 4i + ^15 y si + ^16 y 6i = g li 

y 2i + ^23 y 3i + ^ 24 y 4i + ^25 y 5i + ^26 y 6i " g 2i 

y 3i + ^34 y 4i + ^35 y 5i + ^36 y 6i ~ g 3i 

y 4i + ^45 y 5 i + ^4 6 y 6i “ g 4i 

y 5 i + f 56 y 6i = g 5i 
y 6i = g 6i 
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In the above, the subscript ”i” assumes values from one to seven, corresponding to various so- 
lutions for the seven right-hand side sets of constants. 

The construction of this triangular system of equations is known as the forward solution, 
and the process of obtaining its solution is called back- substitution. The last equation in the 
triangular system gives the value for y 6i directly. If we insert its value in the previous equa- 
tion, we can obtain y s , , and so on for the remainder of the unknowns. Mathematically, the re- 
lation is: 


y 6 i = e 6i 

and 

6 

y j i = i - J2 fikVki 

k = j + 1 

where j = 5, 4, 3, 2, 1 (in that order) and i = 1, 2, . . . , 7 (in any order). 

We have now completed the determination of the Aq = y n by the least-squares principle. 
Theoretically, this procedure may always be followed to a successful conclusion provided that 
the m original equations of condition are independent; that is, provided that the determinant of 
the coefficient matrix C does not vanish. 

The formal solution by the method of least squares is now concluded, but ordinarily a 
measure of the "goodness" of the least- squares fit is desirable. The residuals after the fit 
are assembled in the so-called residual matrix U, equal to B - AX. In terms of elements: 

6 

u i =b. - a., x. (i □ l,2,"* t m) 

i = i 

From this, it is obvious that the sum of the squares (unweighted) of the residuals after the fit 
is given by: 


m 



We can now easily find the so-called variance- covariance matrix of the fit from the in- 
verse C“ 1 of the coefficient matrix in the normal equations. Recall that C“ 1 occupies columns 
two through seven of matrix Y. The variance- covariance matrix is obtained simply by multi- 
plying each element in c~ l by the sum of the squares of the residuals after the fit and dividing 
this product by m-6 (the excess of simultaneous equations of condition over the number of in- 
dependent unknowns). If we represent the variance- covariance matrix by v, then we have that: 


where 


: y. 


i » j + 1 




(i. j = 1,2, 


-. 6 ) 




By comparison with computations performed above in the section titled "The Standard Deviation 
of Fit," we can see that the quantity u { is a standard deviation of fit. More precisely, M f is the 
standard deviation of the residuals after the fit, or the standard deviation of the least squares 
fit. It is not to be confused with 

£ -f ■ /^£ *<“,>’] 
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in the earlier notation, which is the standard deviation of the observational residuals, or the 
standard deviation of the observational fit. 

Finally, we can find the so-called standard errors m, of the six unknowns Aq^y.j . 
These are simply equal to the square-roots of the diagonal elements in the variance-covariance 
matrix, or 


= *V y i.i*i (i = i. 2, .6) 

where . +1 is the term on the principal diagonal of the inverse of the matrix C , corre- 
sponding to the unknown Xj = y n . 

ITERATIVE LEAST-SQUARES PROCEDURE 

The procedure for producing a mean set of orbital elements is essentially an iterated 
least- squares fitting of the first-order Taylor differential expansion of the conditional equa- 
tions to numerous observational values. Using the values for the Aq determined by the method 
of least squares, as described in the preceding section, we can calculate the corrected Izsak 
orbital elements. 

a ' = a + Aq x = a + A a 
e ' - e + Aq 2 = e + Ae 

^0 = ^0 + Ac >3 = ^0 + A7 ?0 

= /3j + Aq 4 = /3j + A/3 t 
fi\ =/3 2 + Aq s =/3 2 +A/3 2 

^3 = ^3 + Ac ^6 = ^3 + A/3 3 

At this point, it is useful to check that the improved or corrected elements are physically mean- 
ingful. For instance, it should be ascertained that the semi- major axis a' > 1 earth radius, that 
the eccentricity e' > 0, that the sine of the inclination 17 ' 0 is not greater than unity in absolute 
value, and so on. 

It is now necessary to update the other parameters used in the differential correction 
process, based upon the improved orbital elements. Accordingly, the various parameters in- 
cluded under the heading, M Prime Constants n” are re-evaluated using the improved set of ele- 
ments. This done, the various parameters included under ’’Mutual Constants” are similarly 
re-evaluated. Now, assuming that the times of the various observations in the data deck are 
available as needed, the Orbit Generator may be used to produce the required calculated values 
of the position and velocity components. From these components, we calculate the local co- 
ordinates of the satellite and then the computed values of the direction cosines (refer to section 
titled, ’’Computation of Direction Cosines”). Finally, the observational residuals are calculated. 
Thus, with each observation time in the data deck are associated corresponding principal time- 
dependent quantities as follows: 

(a) position and velocity components X,Y,Z,x,Y,Z; 

(b) local co-ordinates of the satellite X M ,Y M ,Z M ; 

(c) computed values of the direction cosines L c ,M c ; and 

(d) observational residuals Al,am . 
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Next, the statistical analysis is repeated (refer to the section titled, "The Standard De- 
viation of Fit") wherein the following quantities are determined: the average observational re- 
sidual, the standard deviation of the observational residuals from their mean value, the stand- 
ard deviation of (the observational) fit, the upper and lower range limits for the observational 
residuals, and the standard deviation of fit of the accepted observational residuals. Once these 
quantities are found, the differential correction may be repeated. Of course, the time-inde- 
pendent partial derivatives are computed once only, while the time-varying partial derivatives, 
both with respect to the energy- momenta variables and to the angle-epoch variables, are com- 
puted for each observation time in the data deck. A new set of equations of condition can then 
be assembled and the fitting by least squares repeated. 

In summary, the following sequence of steps represents the iterative least- squares pro- 
cedure in producing a mean set of orbital elements for a given time span represented by a set 
of observation points: 

1. Correct the six Izsak orbital elements utilizing the values determined by the method 
of fitting the equations of condition by least squares. 

2. Update the parameters included in Prime Constants n. 

3. Update the parameters included in Mutual Constants. 

4. Produce sets of position and velocity components for each observation time using the 
Orbit Generator. 

5. Calculate the local co-ordinates of the satellite at each observation time. 

6. Compute the direction cosines of the satellite at each observation time. 

7. Determine the observational residuals for each time point, 

8. Perform a statistical analysis of the observational residuals to find various standard 
deviations and a statistically valid range within which observational residuals must fall for in- 
clusion in the fitting process. 

9. Begin the differential correction process by evaluating the time- independent partial 
derivatives. Then evaluate the time-varying partial derivatives for each observation point. 

10. Assemble the set of equations of condition. 

11. Fit the equations of condition by the method of least squares. First determine the six 
normal equations, then triangularize the system by the Gaussian elimination method, and finally 
use the back- substitution method to find the solution. 

12. Measure the "goodness" of the least- squares fitting by finding the residuals after the 
fit, the variance-covariance matrix, and the standard errors of the unknowns. 

Return to step number one. 


DEFINITIVE ORBITAL PARAMETERS 

The iterative least- squares procedure is generally terminated in one of two ways. Either 
the total number of iterations through the least squares routine is prescribed in advance, or the 
standard deviation of the observational fit is used as the criterion in halting the iterative method. 
If this standard deviation falls below a value prescribed in advance during a given iteration, then 
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the precision of the differential correction is deemed sufficient at that point. Of course, both 
methods of terminating computation can be used concurrently; i.e., if the standard deviation 
does not meet the prescribed criterion by the p-th iteration, then the differential correction 
process is halted. 

At the conclusion of the differential correction, the following definitive orbital param- 
eters may be found: 

The semi-major axis is found by multiplying a by the proper length conversion constant 
(3963.339 mi. /Earth radius or 6378.388 km./Earth radius). 

The eccentricity of the orbit is given by e. 

The inclination of the orbital plane to the equator is given by arc sin v 0 (0° <v 0 < 180°). 

The time of passage through the perigee point is found by multiplying -j3 x by the proper 
time conversion constant (13.4472 min. /Vanguard unit of time or 806.832 sec./Vanguard unit of 
time). The time of perigee passage is given with respect to the reference (or epoch) time t 0) 
which is that used as a basis for the observational times and that corresponding to the initial 
position and velocity components X. ,Y i , z. , X. , Y. , 

The argument of perigee (measured in the orbit plane from the node to the perigee point) 
is found by multiplying by the angular conversion constant 57.295780 deg. /rad. 

The right ascension (measured in the equatorial plane from the vernal equinox) of the 
ascending node is found by multiplying by the angular conversion constant. (Note that these 
last two parameters are angles usually given as greater than or equal to 0° and less than 360°, 
so that some multiple of 360° may have to be added or subtracted to bring the values into this 
principal range). 

The height of the perigee point above the Earth’s surface is found by multiplying a(l-e) - 1 
by one of the length conversion constants given above. 

The height of the apogee point above the Earth’s surface is found by multiplying a(l+e) -1 
by one of the length conversion constants. 

The anomalistic mean motion is found by multiplying a' 3/2 by the angular conversion 
constant and dividing by one of the time conversion constants (this gives the mean motion in 
deg./ min. or deg. /sec.). 

The anomalistic period is found by multiplying 2 77 a 3 ' 2 by one of the time conversion 
constants. 

The mean anomaly (at the time of perigee passage) is found by multiplying a“ 3 2 by 
the angular conversion constant. This expression assumes that the reference (epoch) time t Q 
is zero; in general, the mean anomaly is found by multiplying - a' 3/2 (j3 1 + t 0 ) by the angular 
conversion constant. 


RESULTS OF PRELIMINARY APPLICATIONS 

Both the orbit generator portion and the differential correction process by least-squares 
fitting have been tested independently by application to actual satellite orbits. Primarily, use 
has been made of two relatively close-in yet drag-free satellite orbits, so that neither atmos- 
pheric drag nor luni-solar perturbing forces would exert major disturbing influences. The 
ANNA IB satellite (international designation 1962 BM 1; N.A.S.A. identification number 56017) 
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was launched in October, 1962 under the project direction of the U. S. Navy from the Atlantic 
Missile Range into a near-circular orbit of medium inclination. Its purpose was predominantly 
that of geodetic investigation. The Relay 2 satellite (international designation 1964 3 A; N.A.S.A. 
identification number 64031) was launched in January, 1964 under the project direction of the 
National Aeronautics and Space Administration from the renamed Eastern Test Range into a 
relatively high-eccentricity orbit. Its function was that of active -repeater communications 
satellite. Initial orbital parameters for Doth these satellites are given in Table 1. The obser- 
vational data for the Relay 2 satellite consist of direction- cosine pairs reported from fifteen 
tracking stations in the Minitrack network operated by the N.A.S.A., while the data for ANNA 
IB consist of right ascensions and declinations reported from twelve stations in the optical 
camera network operated by the Smithsonian Astrophysical Observatory. It might be noted that 
no weighting factors were associated with any of the sets of observational data for either the 
Relay 2 or the ANNA IB satellite in the applications described in this section. 


Table 1 

Initial Orbital Parameters for Satellites Used 
in Preliminary Applications 



ANNA IB 

Relay 2 

Perigee (statute miles) 

670 

1298 

Apogee (statute miles) 

728 

4606 

Period (minutes) 

107.8 

194.7 

Inclination to Earth T s 
equator (degrees) 

50.1 

46.0 

Semi-major axis (units of 
Earth’s equatorial radius) 

1.1764 

1.7448 

Eccentricity 

0.00622 

0.23918 

Sine of the inclination 

0.7672 

0.7193 


In order to gauge the intrinsic accuracy of the orbit generator, a double- precision ninth- 
order Cowell numerical integration program was utilized. Two numerically integrated com- 
parison ephemerides were produced: one using recently determined geodetic values for the 
zonal harmonic coefficients in the expansion of the geopotential and the other using the corre- 
sponding values for these coefficients based upon the Vinti potential. Refer to Table 2. The 
numerically integrated ephemeris produced by the geodetic values of the zonal harmonic coef- 
ficients was used as a basis for comparison with both the numerically integrated ephemeris 
produced by Vinti values of the zonal harmonic coefficients and the ephemeris produced by the 
orbit generator based upon the Vinti potential function. Figure 1(a) illustrates the residuals of 
the X-co-ordinate between (1) the Vinti ephemeris and the numerically integrated ephemeris 
produced by geodetic values, and (2) the numerically integrated ephemeris using the Vinti 
values and the numerically integrated ephemeris produced by geodetic values. Figures 1(b) and 
1(c) do likewise for the residuals of the Y-co-ordinate and Z-co-ordinate, respectively. 

The comparisons illustrated in Figure 1 are based on the implicit assumption that the 
initial position and velocity conditions do not contain any inaccuracies. In actual practice, such 
inaccuracies are always present, and they must be removed by utilizing observational data in 
the differential correction. Figure 2 illustrates the determination of a mean set of Izsak orbital 
elements by an iterated least- squares fitting of the differential solution to observational data for 
the ANNA IB satellite. In all cases, the total number of iterations through the least squares 
fitting routine is prescribed in advance to be ten. This number is sufficient to attain conver- 
gence within a very small tolerance. In the graph of each of the six orbital elements, three 
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VINTI EPHEMERIS 


NUMERICALLY 
INTEGRATED 
EPHEMERIS 
USING VINTI 
COEFFICIENTS 


Figure 1 —Co-ordinate residuals for Vinti potential ephemeris and for numerically integrated ephemeris 
using zonal harmonic coefficients of the Vinti potential each compared with numerically integrated 
ephemeris produced by geodetic values. 
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0,7692 
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Figure 2-Convergence of a mean set of Izsak orbital elements in ten iterated least-squares fittings of the differential 
solution to observational data for the ANNA 1 B satellite, for various numbers of equations of condition. 






Table 2 

Zonal Harmonic Coefficients in the Geopotential 
Function Used in Generation of Numerically 
Integrated Comparison Ephemerides 


Coefficient 

Geodetic value 

Vinti potential value* 

J 2 

1.0823 x 10- 3 

1.0823 x 10- 3 

J 3 

-2.3 x 10- 6 

0 


-1.8 x 10- 6 

-1.2X10- 6 

J (n > 5) 

<1 x 10- 6 

< 1 X 10- 8 " 


♦The zonal harmonic coefficients for the Vinti potential function are obtained from the 
relations: J 2k = (-l) k+1 J k and J 2k+1 = 0. 


"curves" (actually a sequence of connected line segments) are shown, corresponding to various 
numbers of observations included in the fitting. An equation of condition results, of course, 
from a "semi-observation": either a single right ascension or a single declination value. One 
curve represents the minimum number of equations of condition for a true least- square fitting, 
viz., seven. This is associated with an observational arc length of approximately 45 hours. A 
second curve represents twenty equations of condition, or an addition of eleven equations, ex- 
tending the observational arc length to approximately 75 hours. The third curve represents 
fifty equations of condition, or a further addition of thirty equations, extending the observational 
arc length to a total of approximately 98 hours. The starting point of each of the three arcs is 
the same, so that they overlap in time. Notice that each observational arc produces a somewhat 
different set of mean orbital elements, depending upon the additional observational values intro- 
duced. Physically, this may be explained as the resultant effect of forces not accounted for in 
the analytical theory. For example, electromagnetic disturbances, solar radiation pressures, 
aerodynamic drag, meteoric bombardment, etc., all influence the mean set of orbital elements 
to the extent that they are reflected in the observational values. In performing the iterated 
least-square fittings, all the residuals corresponding to the pre-selected observation times 
were accepted at each fitting. That is, the acceptable range of values for the residuals con- 
stituted infinitely wide bands on either side of the mean value of the residuals. Mathematically, 
using symbols introduced in the section titled, "The Standard Deviation of Fit," 

t r r r 2 l “ 1 im [R - jo*. R + j cr] . 

j-00 

Figure 3 illustrates the determination of a mean set of Izsak orbital elements by an iter- 
ated least- squares fitting of the differential solution to observational data for the Relay 2 satel- 
lite. In this case, however, the observational arc length and the total number of observations 
are held fixed, while the acceptability criterion for the observational residuals is varied. The 
arc length in all cases is one week, representing a total of eighty observations or a maximum 
of 160 possible conditional equations. In each of the six graphs, one curve corresponds to a 
"three- sigma" criterion, i.e., 


[ r r r 2 ] - [R - 3cr, R + 3cr] 

where, of course, cr is the standard deviation of the observational residuals from their mean 
value. A second curve corresponds to a two-sigma criterion, and the third curve to a one- 
sigma criterion. Each curve is terminated when convergence of the orbital element is attained. 
Notice that convergence appears to be a slower process with a one- sigma criterion than with 
either a two- sigma or a three- sigma criterion. The rate of convergence in these latter two 
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Figure 3-Convergence of a mean set of Izsak orbital elements in ten iterated I east -squares fittings of the 
differential solution to observational data for the Relay 2 satellite, for various observational residual 
acceptability half-width criteria. 





cases seems generally about the same. One effect of a wider acceptance range appears to be 
greater fluctuations in the value of an orbital element early in the iterated fitting procedure, 
although this is not always true. Also, despite the fact that differing numbers of conditional 
equations are accepted in the fittings depending upon the criterion for the residuals, the values 
of the orbital elements at convergence are remarkably similar. Refer to Table 3 for precise 
values, including the required number of iterations to attain convergence in each case. The 
uncertainties in the final significant figures (stated as ±X) are estimates based upon slight 
fluctuations in the values of the orbital elements in least- square fittings after convergence is 
attained. 


Table 3 

Values at Convergence of Izsak Orbital Elements for Varying 
Observational Residual Acceptability Half-Widths 


Orbital Elements 

One Sigma 
Criterion 

Two Sigma 
Criterion 

Three Sigma 
Criterion 

Semi-major axis a 

1.7444277 (12) 

1.7444278 
± 0 

1.7444278 , . 
± 0 yi) 

Eccentricity e 

0.2391624 , 01 , 
± 4 U1 ' 

0.2391728 . . 
± 3 (,} 

0.2391818 . . 
± 2 (l> 

Sine of inclination rj Q 

0.7231110 

± 2 ’ 

0.7231020 , . 
± 4 

0.7231015 ... 
± 2 w 

Time of perigee passage p x 

0.099885 , 00 . 

± 6 (ZZ) 

0.099985 /qv 
± 7 W 

0.099942 
± 9 (8) 

Argument of perigee p 2 

3.222265 , 01 . 

± 5 

3.222232 ... 

± 5 w 

3.222278 , . 

± 6 

Rt. asc. of ascending node fi 

-2.380274 
± 3 

-2.380244 (8) 

± 4 ' ' 

-2.380241 ... 

± 4 (8) 


Note: Units of all elements are canonical (a in Earth equatorial radii; p j in Vanguard units of time; ^ an ^ °3 radians). 
The integers in parentheses refer to the number of iterations required to attain the converged value given. 


Table 4 presents the same information relative to the orbital elements as Table 3, but 
for an acceptance criterion fixed at two sigma, with the maximum possible number of condi- 
tional equations varied. The observational arc length remains one week, but the 160 maximum 
possible number of conditional equations are first reduced to one hundred, and then this num- 
ber is in turn reduced to forty. An attempt was made to maintain an even distribution of the 
observations throughout the seven- day period, while still operating on the " subset principle" 
(i.e., the set of twenty observations is a subset of the set of fifty observations, which is in turn 
a subset of the original set of eighty observations). 

Table 5 again records the same information relating to the orbital elements, but this time 
the parameter involving the order of precision in the differential correction is varied. Here the 
maximum possible number of conditional equations covering the one-week observational arc is 
held constant at forty, and the acceptance criterion is fixed at two sigma. The inexact designa- 
tions "first order" and "second order" indicate whether or not terms of purely second order 
are retained in the differential correction. (Refer to the section titled, "Analytical Procedure 
of Differential Correction.") It is seen that retaining terms of purely second order adds im- 
measurably to the precision of the final converged results in all cases, and, similarly, does not 
affect the rate of least-squares convergence. 
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Table 4 

Values at Convergence of Izsak Orbital Elements for Varying 
Numbers of Observational Points 


Orbital Elements 

40 Conditional 
Equations 

100 Conditional 
Equations 

160 Conditional 
Equations 

Semi- major axis a 

1.7444276 , 1rt , 

± o 110) 

1.7444279 
± 0 ' lu; 

1.7444278 
± 0 (i) 

Eccentricity e 

0.2391425 . . 

+ 5 

0.2391736 /.v 
± 2 ' y ' 

0.2391728 . . 
± 3 ''' 

Sine of inclination rj Q 

0.7231009 , 10 . 
± 3 (1Z) 

0.7230975 . m 
± 3 UU; 

0.7231020 . . 
± 4 { ’ 

Time of perigee passage p x 

0.100.16 (12) 

0.100002 
± 7 ui; 

0.099985 , Q > 
± 7 w 

Argument of perigee ;3 2 

3.22 2 097 , u 
± 5 

3.222261 . 0 . 

± 6 w 

3.222232 . . 
± 5 ( ’ 

Rt. asc. of ascending node p 3 

-2.380249 (12) 

-2.380248 . . 

± 3 UU; 

-2.380244 
± 4 


Note: Units of all elements are canonical (a in Earth equatorial radii; p x in Vanguard units of time; fl 2 and in radians). 
The integers in parentheses refer to the number of iterations required to attain the converged value given. 


Table 5 

Values at Convergence of Izsak Orbital Elements for Varying 
Orders of Precision in Differential Correction 


Orbital Elements 

First Order 

Second Order 

Semi-major axis a 

1.7444276 {1Q) 

1.7444276 nm 

± o u ' 

Eccentricity e 

0.2391427 . . 

± 8 U1; 

0.2391425 (n) 

Sine of inclination r / 0 

0.7231010 
± 3 

0.7231009 
± 3 {l£) 

Time of perigee passage ? x 

oaoons (12) 


Argument of perigee 2 

3.222098 . . 

± 6 

3.222097 . . 

x 5 (11) 

Rt.asc 0 of ascending node - 3 

-2.380250 (12) 

-2.380249 no . 
x 3 [ > 


Note: Units of all elements are canonical (a in Earth equatorial radii; pj in Vanguard units of time; v, 
and ^ in radians). The integers in parentheses refer to the number of iterations required to attain 
the converged value given. 
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The remainder of the Figures display the convergence of perhaps the most significant 
single parameter in evaluating the efficacy of the differential correction process, viz., the 
standard deviation of fit. Actually, there are two standard deviations shown in each graph. 

The upper curve corresponds to a standard deviation of fit which includes all of the observa- 
tional residuals, while the lower curve corresponds to a standard deviation of fit which in- 
cludes only the observational residuals accepted at each fitting. Plotted on the same abscissa 
is a curve showing the number of equations of condition (or, equivalently, the number of ob- 
servational residuals) accepted at each iteration of the fitting process. 

Figure 4 illustrates the standard deviations for a maximum of forty possible conditional 
equations covering a one-week observational arc for the Relay 2 satellite. Note that conver- 
gence using a two- sigma criterion for the residuals, as shown in Figure 4(b), is much more 
rapid than the convergence using a one-sigma criterion shown in Figure 4(a). However, the 
convergence is not so smoothly monotonic in the case of the wider acceptance range. Both 
these facts confirm what was said earlier about the convergence of the orbital elements. 

Figure 5 illustrates the standard deviations for a maximum of one hundred possible con- 
ditional equations and Figure 6 for a maximum of 160 possible conditional equations, both cov- 
ering the same one- week observational arc for Relay 2. Similar remarks apply to these Fig- 
ures as to Figure 4. Table 6 supplies the values of the standard deviations at convergence for 
the various runs illustrated in Figures 4, 5, and 6, as well as several others not graphed. It 
also gives the number of accepted residuals at convergence, and the number of iterated fittings 
required to achieve convergence in each case. 

Table 6 

Values of Standard Deviations of Fit and Number of Accepted 
Conditional Equations at Convergence for Various Runs 
Covering a One-Week Observational Arc 


Description of Run 

Std. Dev. 
of Fit 
(all) 

Std. Dev. 

of Fit 
(accepted) 

Accepted 

Residuals 

Percentage 
of Total 

Iter- 

ations 

Total 

Conditional 

Equations 

Residual 

Criterion 

Order 
of D.C. 

40 

ler 

2nd 

0.432 

0.157 

37 

92.5 

26 

40 

2 a 

2nd 

0.440 

0.177 

38 

95 

12 

40 

2a 

1st 



38 

95 

12 


la 

2nd 



95 

95 

11 


2a 

2nd 

0.377 

0.169 

97 

97 

9 


la 

1st 

0.307 


141 

88.1 

24 

160 

la 

2nd 

0.307 


140 

87.5 

22 

160 

2 a 

2nd 

0.313 

0.172 

157 

98.1 

8 

160 

3a 

2nd 

0.315 

0.187 

158 

98.75 

9 


Note: All standard deviations of fit are given in mils (i.e., in units of 10~ 3 ). The parenthetical word "all” signifies that 

all of the observational residuals, A L and AM, were included in determining the standard deviation of fit; "accepted” 
means that only the observational residuals corresponding to the accepted conditional equations were included in 
determining the standard deviation of fit. 
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Figure 6— Standard deviations of the observational residuals and the number of equations of condition 
accepted at each iteration of the fitting process for a maximum of 160 possible conditional equations 
covering a one-week observational arc for the Relay 2 satellite. 
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Figure 7 illustrates the standard deviations for a maximum of one hundred possible con- 
ditional equations covering an observational arc of only three hours for Relay 2. This is the 
three-hour period immediately following insertion of the satellite into orbit, when observations 
are recorded at very frequent intervals in order to insure a wealth of data for the real-time 
differential correction. Here, using a one-sigma criterion, convergence of the orbital elements 
occurs after only four (in some cases, five) iterations. The standard deviations of fit converge 
after three iterations to values of 0.425 x 10~ 3 (all one hundred observational residuals) and 
0.145 x 10 " 3 (including seventy-seven accepted observational residuals). The graph shows that 
the standard deviations remain essentially constant after the third iteration, and this is con- 
firmed by the insignificant fluctuations in the orbital elements after the third iteration, although 
a total of ten iterations through the least squares fitting routine was prescribed in advance. 

Figure 8 illustrates the standard deviations for a maximum of one hundred possible con- 
ditional equations for two distinct non-overlapping observational arcs for the ANNA IB satel- 
lite. Figure 8(a) covers an arc of approximately nine days and fifteen hours, while Figure 8(b) 
covers an arc of approximately six days and nineteen hours. Both use a one-sigma criterion, 
and convergence of the standard deviations occurs after five iterations in both cases. The 
values are a relatively large 38.379 milliradians (all one hundred observational residuals) and 
10.919 mrad (including eighty-eight accepted observational residuals) for Figure 8(a). For the 
somewhat shorter arc in Figure 8(b), the values are 6.684 mrad (all one hundred observational 
residuals) and 0.726 mrad (including ninety-five accepted observational residuals). 

t 

The totality of data presented herein represents a small sampling of the preliminary 
applications by which the orbit generator and differential correction have been tested. Yet 
this sampling is indicative of the utility of the spheroidal method for artificial satellite orbits. 


CONCLUDING REMARKS 

The method of solution for unretarded satellite orbits discussed in this paper has been 
programmed, primarily in the FORTRAN language, for use on the I.B.M. 7094 digital electronic 
computer. It requires a relatively small number of computer storage locations, and the ana- 
lytical nature of the entire procedure assures a very rapid computational process. Extensive 
tests have indicated a capacity for generating co-ordinate and velocity points, based upon a set 
of empirically estimated initial conditions, in impressively short intervals of computer oper- 
ating time. 

Presently, work is underway on slightly modifying the accurate reference orbit to account 
for the effects of the most important perturbations of the neglected zonal harmonics, notably 
the third and the residual fourth. The inclusion of these perturbative effects by a procedure 
described in a recent paper by Vinti (see references) is expected to improve the accuracy of 
the method so as to provide computed values agreeing with observation over a longer interval 
of time. 

In the future, a method of modifying the spheroidal potential for an oblate planet in order 
to permit the exact inclusion of the effects of the third zonal harmonic in the reference orbit is 
anticipated. Preliminary investigations are also being conducted into accounting for the luni- 
solar forces and aerodynamic drag. Further results will be published as they become available 
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Figure 7 — Standard deviations of the observational residuals and the number of equations of 
condition accepted at each iteration of the fitting process for a maximum of 100 possible 
conditional equations covering a three-hour observational arc for the Relay 2 satellite. 
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Figure 8-Standard deviations of the observational residuals and the number of equations of condition 
accepted at each iteration of the fitting process for a maximum of 100 possible conditional equations 
covering two distinct observational arcs for the ANNA 1 B satellite. 
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APPENDIX 


Herein we present the modifications which must be introduced in order to utilize an 
alternate form of satellite tracking data known as right ascension-declination data. Such data 
are recorded, for instance, by the optical Baker-Nunn cameras of the Astrophysical Observ- 
atory of the Smithsonian Institution. The modifications to be described replace the material 
presented in the main body of this report in the section titled, "Computation of Direction 
Cosines." 

A set of observation data of the right ascension-declination type includes the following 
parameters for each recorded spacecraft observation: 

t' , the date and time of observation. The same remarks about removing reference to 
the calendar in transforming t' to the relative time t apply here as included in the main body 
of this report. 

k, the code number for the tracking station reporting the observation. 

a 0 , the observed right ascension, measured in radians eastward from the vernal equinox 

(0 < a 0 < 2 77 ). 

S 0 , the observed declination in radians, measured as positive north of the equator and 
as negative south of the equator (-tt/2 < 8 0 < + n/ 2). 

w a and w s , the weighting factors corresponding to observations a Q and s 0 , respectively. 
This information is optional; if not provided, then it is assumed that w a and w s are each unity. 

The co-ordinate system employed for the observation data is centered at the tracking 
station on the Earth 1 s surface, and, unlike the system used for recording direction- cosine 
data, its three co-ordinate axes are parallel to the respective axes of the inertial system. 

That is, the Z-axis is parallel to the Earth’s polar axis, and the X-Y plane is parallel to the 
equatorial plane of the Earth, with the X-axis extending toward the vernal equinox. The Y- 
axis extends orthogonally to the east to form a right-handed system. 

The differential correction process requires the same data to be available as listed in 
the main body of this report, viz., the Earth’s flattening coefficient f, the Earth’s rotational 
rate a>, the geodetic longitudes \ E of the stations, the geodetic latitudes 6 h of the stations, the 
altitudes H of the stations, the angular distances from the vernal equinox to the Greenwich 
meridian at midnight Greenwich mean time for each day in the observational arc, and the 
reference time t 0 . 

Computations follow the same scheme given in the main body of this report for the follow- 
ing parameters: the auxiliary functions c and S, the geocentric latitude 0 Q , the geocentric 
distance p of the station, the angular distance 8 between the vernal equinox and the observa- 
tion meridian plane, and the inertial geocentric co-ordinates x T , Y T , Z T of the station. How- 
ever, the angle \p x between the vernal equinox and the tracking station’s X-co-ordinate axis 
(measured in the observation latitude plane) is zero, since the topocentric and inertial co- 
ordinate systems are parallel. No rotations are necessary to bring the two systems into co- 
incidence; a single translation will suffice. Hence, the relations for the topocentric or local 
co-ordinates of the satellite are simply: 
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X* =X-Xr 


Y . = Y - Y t 
\ = Z - Zt 

where X, Y, Z are the inertial geocentric co-ordinates of the satellite predicted by the orbit 
generator. The above simplified relations are obtained from those of the direction-cosine- 
data case by the artificial device of setting ^ = 0 and 0 B = n/2 in the corresponding equa- 
tions for X M , Y„, Z M given in the main body of this report. (Refer to the note at the end of this 
appendix.) 

The computed values of the right ascension and the declination may now be found in terms 
of the local co-ordinates: 


a - 
c 



arctan 


(Xj; + Yj*) I/2 


It is important that the angles a c and b c be placed in the proper quadrant for comparison 
purposes with the angles a Q and 8 0 . In the case of the right ascension, this is done by examin- 
ing the signs of and y m separately. The following list presents all possible combinations 
(note that the range for a Q is 0 < a 0 < 2 77). 


X u > 0, Y u > 0 : 0 < a < — 

•I’M co 


X M > 0, Y m < 0 : ~ < a c < 2v 


X u > °, Y m □ 0 : 


a - 0 


X M < °, Y m > 0 : 


< a <77 


X M < 0, Y u < 0 : n < a < — 


X„ < 0, Y M = 0 : 


0, Y„ > 0 : 


0, Y„ < 0 : 


2 

377 


0, Y u = 0 : a indeterminate 

M C 


In the case of the declination, the signs of the numerator and denominator of the arctangent 
argument are examined separately. The following list presents all possible combinations (note 
that the range for S Q is - tt/2 < < + tt/2). 
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( x m + YJ ) 1/2 > 0 , 

Z M > 0 : 

0 < s < — 

c 2 

+ Y S > 1/2 > o . 

Z M < 0 : 

--<8 <0 

2 

( X M + Y m ) 1/2 > 

II 

o 

8 c = 0 

( x « + y m > ,/2 = 0 . 

Z M > 0 : 

8 - — 

c 2 

( x £ + y m ) 1/2 = 0 , 

Z M < 0 : 

8 = -~ 
c 2 

,) 1/2 = Z M - 0 is not physically possible 


The observational residuals are now found: 


Aa - a n - a 

0 c 


AS = 8. - 8 

0 c 


Here too, care must be exercised. There is one instance where simple subtraction in finding 
the observational residual will yield a misleading result. If one of the right ascensions 
(either observed or computed) is in the first quadrant and very nearly zero and the other right 
ascension is in the fourth quadrant and very nearly 2 rr , then direct subtraction will provide an 
erroneous result near to 2 77, whereas the intended difference is near to zero. This situation 
can be rectified by the following logical steps: 

If |a 0 - aj < 77, then Aa = a 0 - a c (as above). 

If I a Q - a c | >77, then Aa = Sgn (a 0 - a c ) [277 - | a Q - a c | ] . 

Equivalently, whenever |a 0 - aj > 77 , use the following: 

(1) if a > a , then Aa - 277 - a n + a >0. 

u c U c 

(2) if a < a , then Aa = a _ a. - 277 < o. 

The statistical analysis of the observational residuals follows the procedure given in the 
main body of this report in the section titled, ’’The Standard Deviation of Fit” except that the 
observational residuals are given by Aa. and as., rather than by AL a and AM t . Hence, the 
average residual is given by: 


n 



The standard deviation of the residuals from their mean value is found from: 


< 7 



2^ [(Aa. -R) 2 + (AS. -R) 2 ] 

i=l 
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The standard deviation of fit is given by: 


cr 


f 




Modifications are now presented to supplant the material from the main body of this re- 
port in the section titled, "Analytical Procedure of Differential Correction." 

The first-order Taylor series expansion of the equations of condition may be written: 


Act 


f-’K 

4 r 3q - 


A< ii 


AS 


V s -! 38 

i = l 1 


where q. (i = 1, 2, ••• , 6) are the mean or Izsak orbital elements. Expanding the above 
partial derivatives by the chain rule yields: 

Ba c da c BXjj Ba c BY M Ba c 

3X„ 3q. + 3Y„ dq. + 3Z„ 3q. 

db c ?S c dX u dS c 3Y„ B8 c 3Z M 

” ?X M + BY M H + ^ 

From the equations for a c and S c in terms of the local co-ordinates, we find: 


da 


Y « <XS - Y 2 )' 1 



+ \ (% + Y 2> 


da 




38 


- X m Z m (X 2 + Y^)- 1/2 (X* 


Y 2 


z«)' 
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B 5 

3^ = - Y « Z « <*S + Y^)' 1/2 (X* + Y* + Z2)- 1 

It = + (34 + Y^) I/2 ( 3 ^ + Y * + Z*)- 1 


Since the station co-ordinates X T , y t , z t are independent of orbital parameters (and merely 
geodesic functions), the following simple relations hold: 

^ Y m _ _BY_ 

BZ m _ BZ 
^ 3 q . 


The method for calculating the partial derivatives Bx/Bq. , By/Bq., and B^/Bq. is identical to 
that presented in the differential correction scheme in the main body of this report* Then the 
equations of condition are formulated in a precisely analogous manner to that given for the 
direction- cosine data (see the section titled, "The Equations of Condition"), and there is little 
need to repeat the explicit form of these equations* 

NOTE: The fact that the observational co-ordinate system is independent of the latitude 

and longitude of the tracking station for right ascension-declination data (as is not the case for 
direction-cosine data) leads to certain possible simplifications in the determination of the 
computed co-ordinates of the satellite, a c and S c * First, recall that the equations for the 
Cartesian inertial co-ordinates of the observation point are given by: 

X T = p cos 9 q cos 8 
- P cos 0 Q sin S 
Z T = p sin 9 g 

where 5 = (\ 0 ) d + w(at) + \ E . Here the terms (\ 0 ) d and a>( AT) depend upon the time of the 
observation only, while the term \ E is a function of the location of the observation point* Let 
us denote: 


8' = (\>) d + *>(AT) 


Then we can expand the above equations as: 

X T = p cos 6 q cos (S' + \ £ ) 

= p cos cos )v E cos S' — p cos Oq sin A. £ sin 8 1 
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Y T - P cos 6 q sin (S' + A £ ) 

- P cos & c cos Aj sin S' + p cos sin cos S' 




Now denote : 


so that: 


Z T ~ p sin d Q 


Xq = p cos 6 q cos Aj. 

Y q = p cos 0 Q sin 
Z Q □ p sin 0 Q 



cos S' -sin S' 0 


X 0 

sin S' cos S' 0 


Y o 

0 0 1 


Z °_ 


This represents, in matrix form, the fact that the co-ordinates x T , Y T , z T are obtained from 

x 0 , Y 0 , z 0 by a simple rotation about the inertial Z-axis through an angle S'. Here s' is the 

angle between the vernal equinox and the Greenwich meridian at observation time. The rec- 
tangular co-ordinates x 0 , Y 0 , z 0 , obtained directly from the spherical geocentric co-ordinates 
P, e c > station, represent the Cartesian inertial geocentric co-ordinates of the track- 

ing station at a time when the Greenwich meridian and the first point of Aries (the vernal 
equinox) coincide. If the co-ordinates X 0 , Y 0 , z 0 (dependent upon the station location only) are 
provided as input parameters rather than a £ , 0 d , and h, then the computations leading up to 
X T , y t , Z t are simplified considerably. We need not first compute C, s, 6 Q9 p, and S. In- 
stead, find s' from parameters relating to the time of observation, and then compute directly: 

X T = X 0 cos S' - Y 0 sin S' 

Y t = X 0 sin S' + Y 0 cos S' 

Zt = Z 0 

Note that this simplified procedure cannot be adopted efficiently with direction-cosine data 
because the rotation matrix involved in computing X T , Y T , z T is a function of # D , the geodetic 
latitude, and 0 x , an angular parameter dependent upon the geodetic longitude. 
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